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1  Summary  of  the  Research  Results 


The  following  goals  have  been  accomplished  in  the  course  of  this  research  project:  A  method 
for  large  eddy  simulations  which  estimates  subgrid  scales  in  terms  of  the  resolved  scales  was 
developed.  In  the  procedure  the  resolved  large  scale  velocity  field  is  expanded  to  subgrid 
scales  two  times  smaller  than  the  LES  grid  scale  by  utilizing  properties  of  a  filtering  procedure 
and  the  representation  of  quantities  in  terms  of  basis  functions  such  as  Fourier  or  Legendre 
polynomials.  The  phases  associated  with  the  newly  computed  smaller  scales  axe  then  ad¬ 
justed  in  order  to  correspond  to  the  small  scale  phases  generated  by  nonlinear  interactions 
of  the  large  scale  field.  The  subgrid-scale  term  in  the  LES  equation  can  then  be  evaluated 
directly  using  the  modified  velocity  field.  This  approach  to  determine  the  subgrid-scale  term 
naturally  accounts  for  backscatter  and  does  not  require  any  adjustable  constants. 

The  modeling  procedure  was  implemented  for  isotropic  turbulence  and  turbulent  channel 
flow.  In  a  priori  tests  for  isotropic  turbulence  databases  correlation  coefficients  between 
modeled  and  exact  stresses  exceeded  90%  and  the  modeled  mean  and  rms  values  were  within 
5%  of  the  exact  values  (for  nonvanishing  components).  The  level  of  agreement  between 
the  modeled  and  the  exact  SGS  quantities  is  significantly  better  than  for  typical  models 
currently  in  use,  and  is  achieved  without  the  presence  of  any  adjustable  constants.  Forward 
and  backscatter  are  also  predicted  in  a  very  good  agreement  with  the  exact  values.  A  priori 
analysis  of  the  channel  flow  databases  for  Rer  =  180  gave  equally  encouraging  results. 

The  model  has  also  been  implemented  in  the  actual  large  eddy  simulations  of  channel 
flow  at  two  different  Reynolds  number  and  was  compared  with  the  Smagorinsky  model 
optimized  for  this  case  as  well  as  with  the  exact  values  computed  from  the  well  resolved 
DNS.  For  all  quantities  the  prediction  of  the  new  model  is  much  better  than  obtained  with 
the  Smagorinsky  model.  In  particular  plane  averages  of  the  diagonal  components  of  the  SGS 
stress  tensor  for  the  Smagorinsky  model  are  essentially  equal  to  zero  while  the  DNS  results 
and  the  new  model  give  non-zero  results  in  good  mutual  agreement.  The  only  non-zero 
component  ri3  for  the  Smagorinsky  model  is  significantly  underpredicted  compared  with 
the  DNS  data  and  the  new  model.  Very  similar  results  would  be  expected  if  the  dynamic 
model  were  used  in  place  of  the  Smagorinsky  model  because  of  the  basic  assumption  in  both 
models  of  the  proportionality  between  stresses  and  rates-of-strain.  In  all  investigated  cases 
the  conclusion  is  that  the  estimation  model  is  very  robust  and  provides  results  in  excellent 
agreement  with  the  exact  numerical  or  experimental  data  and  generally  performs  much  better 
than  typical  SGS  models  currently  in  use.  Among  its  strengths  is  lack  of  adjustable  constants, 
presence  of  the  local  backscatter  which  does  not  cause  any  numerical  instabilities  so  that  ad 
hoc  averaging  procedures  can  be  avoided,  and  by  construction  the  Galilean  invariance  and 
the  proper  near-wall  behavior. 
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4  Research  Results 


This  section  consists  of  submitted  journal  papers  which  provide  detailed  information  about 
research  results  summarized  in  section  1. 
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A  subgrid-scale  model  based  on  the  estimation  of  unresolved  scales  of  turbulence 


J.  Andrzej  Domaradzki  and  Eileen  M.  Saiki 

Department  of  Aerospace  Engineering, 

University  of  Southern  California, 

Los  Angeles,  CA  90089-^1191,  U,S,A,, 
tel  213-740-5357,  FAX  213-740-7774 
(October  15,  1996) 

A  new  method  for  large  eddy  simulations  is  described  and  evaluated.  In  the  proposed  method  the 
primary  modeled  quantity  is  the  unfiltered  velocity  field  appearing  in  the  definition  of  the  subgrid- 
scale  stress  tensor.  An  estimate  of  the  unfiltered  velocity  is  obtained  by  exp€mding  the  resolved  Ifirge 
scale  velocity  field  to  subgrid-scales  two  times  smaller  than  the  grid  scale.  The  estimation  procedure 
consists  of  two  steps.  The  first  step  utilizes  properties  of  a  filtering  operation  and  the  representation 
of  quantities  in  terms  of  basis  functions  such  as  Fourier  polynomials.  In  the  second  step,  the  phases 
associated  with  the  newly  computed  smaller  scales  are  adjusted  in  order  to  correspond  to  the  small 
scale  phases  generated  by  nonlinear  interactions  of  the  large  scale  field.  The  estimated  velocity  field 
is  expressed  entirely  in  terms  of  the  known,  resolved  velocity  field  without  any  adjustable  constants. 
The  modeling  procedure  is  evaluated  in  a  priori  analyses  using  direct  numerical  simulation  results 
of  channel  flow  at  low  Reynolds  number  and  in  actual  large  eddy  simulations  of  channel  flow  at 
two  different  Reynolds  numbers.  In  all  cases,  the  new  model  performs  better  than  or  comparable 
to  classical  eddy  viscosity  models  for  the  majority  of  physical  quantities.  In  particular,  the  subgrid- 
sccde  stress  tensor  is  predicted  very  accurately  and  the  procedure  naturally  £iccounts  for  backscatter 
without  any  adverse  effects  on  the  numerical  stability. 

47.27.Eq 
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1.  INTRODUCTION 


The  large  eddy  simulation  equations  for  an  incompressible  flow  are  conventionally  written  as 

d_  d _  I  d  _  d 


_8_ 

dxi 


Ui  =  0. 


The  overbar  denotes  spatial  filtering  which,  for  a  quantity  /(x),  is  defined  by  the  integral 


(2) 


/(x)  =  J  G(x,x')/(x'),  (3) 

where  G  is  a  given  kernel  function.  Using  the  filtering  operation  (3),  /  is  decomposed  into  a  large  scale  (or  resolved) 
component  /  and  a  subgrid-scale  (or  unresolved)  component 

/(x)  =7(x)  +  /'(x).  (4) 


In  Eq.  (1),  Ui,  p,  and  1/  are  the  velocity,  pressure,  and  the  kinematic  viscosity,  respectively,  and  the  effects  of 
subgrid-scale  quantities  on  the  resolved  velocity  are  described  by  a  subgrid-scale  stress  tensor 


Tij  =  UiUj  -  Ui  Uj,  (5) 

which  must  be  modeled  in  terms  of  the  resolved  quantities  to  close  Eq.  (1). 

Currently,  several  subgrid-scale  (SGS)  models  are  commonly  used  in  large  eddy  simulations  (LES)  of  turbulence. 
Recent  reviews  of  SGS  modeling  are  given  by  Piomelli  [1]  and  Lesieur  [2]  and  many  applications  are  described  in  the 
Proceedings  edited  by  Galperin  and  Orszag  [3].  The  SGS  models  fall  into  three  general  categories:  eddy  viscosity 
models,  similarity  models,  and  so-called  mixed  models  which  combine  eddy  viscosity  and  similarity  expressions. 

The  eddy  viscosity  models  assume  that 


Tij  =  2UtSij  +  ^TkkSij, 


where  Sij  is  the  resolved  rate-of-strain  tensor 


S 


ij 


1 

2 


fduj  duj\ 
\dXj^  dXi)^ 


(6) 

(7) 


and  Ut  is  the  eddy  viscosity.  The  trace  term  in  (6)  is  always  combined  with  the  pressure  term  in  Eq.  (1).  Among 
the  eddy  viscosity  models,  the  Smagorinsky  model  [4]  and  the  models  for  homogeneous  turbulence  proposed  by 
Kraichnan  [5]  and  Chollet  and  Lesieur  [6]  have  a  status  of  classical  models.  These  models  are  also  used  as  a  starting 
point  to  construct  more  recent  models:  the  dynamic  model  of  Germano  et  al.  [7]  and  its  variations  [8-10],  and  the 
structure  function  model  of  Metais  and  Lesieur  [11].  The  most  common  and  important  feature  of  the  SGS  eddy 
viscosity  models  is  that  they  properly  account  for  global  subgrid-scale  dissipation,  i.e.  the  net  energy  flux  from  the 
resolved  to  the  unresolved  subgrid-scales.  This  feature  results  in  good  predictions  of  important  turbulent  quantities 
such  as  mean  velocities  and  rms  velocity  fluctuations.  On  the  other  hand,  the  models  often  fail  in  predicting  other 
quantities.  For  example,  in  a  priori  analyses  of  direct  numerical  simulations  (DNS),  correlation  coefficients  between 
the  actual  and  modeled  SGS  quantities  (stresses,  dissipation,  nonlinear  term,  etc.)  are  very  low,  usually  between  0.1 
and  0.3  [12-14].  Moreover,  the  eddy  viscosity  formulation  is  inherently  ill-suited  to  model  the  phenomenon  of  inverse 
energy  transfer  from  the  subgrid-scales  to  the  resolved  scales.  This  process,  frequently  referred  to  as  backscatter, 
may  be  quite  significant.  In  analyses  of  direct  numerical  simulations  of  various  turbulent  flows  [14-17],  backscatter 
was  found  to  be  comparable  to  or  often  larger  than  the  net  subgrid-scale  dissipation.  The  Smagorinsky  and  the 
structure  function  models  are  unable  to  predict  backscatter,  and  attempts  to  model  backscatter  in  the  framework 
of  the  dynamic  model  may  lead  to  numerical  instabilities,  because  the  dissipative  Smagorinsky  expression  is  used  to 
describe  non-dissipative  backscatter  effects.  Additional  procedures  limiting  excessive  backscatter  are  needed  to  prevent 
the  numerical  instabilities  in  this  latter  case  [10,18].  In  order  to  model  backscatter  effectively,  the  common  practice 
is  to  introduce  an  additional  forcing  term,  which  is  either  random  or  deterministic,  in  the  LES  equations  [19-22]). 
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Backscatter  is  naturally  present  in  similarity  models.  Similarity  models  [23,24]  assume  that  the  unknown  subgrid- 
scale  stress  tensor  can  be  approximated  by  a  stress  tensor  calculated  from  the  resolved  field  employing  an  additional 
filtering  with  the  filter  width  equal  to  or  larger  than  the  one  used  to  obtain  the  original  resolved  field.  The  Bardina 
model  is  given  by  the  following  expression 


TijfauiUj-UiUj.  (8) 

The  similarity  model  of  Liu  et  al,  [24]  assumes  that  the  SGS  stress  tensor  is  proportional  to  the  resolved  stress  tensor 
Lij  which  is  computed  from  the  resolved  velocity  using  a  wider  filter  (denoted  by  hat) 

(9) 

The  stresses  calculated  in  the  similarity  models  are  found  to  correlate  very  well  with  the  exact  stresses  in  a  priori 
analyses.  However,  because  of  the  additional  filtering  the  modeled  SGS  interactions  affect  scales  larger  than  those 
affected  by  the  exact  SGS  interactions.  This  feature  appears  as  a  mismatch  in  the  characteristic  length  scales  for  the 
modeled  and  the  exact  SGS  fields  [24].  Also,  the  model  of  Bardina  [23]  significantly  under  predicts  net  SGS  dissipation 
and  consequently,  it  cannot  be  used  to  reliably  predict  mean  and  rms  quantities  in  actual  large  eddy  simulations.  In 
particular,  the  SGS  stress  tensor  in  the  Bardina  model  vanishes  if  the  sharp  spectral  filter  is  employed.  In  a  priori 
analyses,  the  similarity  model  proposed  by  Liu  et  al,  [24]  provides  correct  amounts  of  SGS  dissipation  in  addition  to 
good  correlations  between  the  modeled  and  exeu:t  stresses.  However,  when  the  model  is  used  in  actual  LES  it  appears 
to  experience  the  same  difficulties  as  the  Bardina  model  (C.  Meneveau,  private  communication). 

Mixed  models  [23,25,26]  attempt  to  combine  the  good  dissipative  features  of  the  eddy  viscosity  models  with  the 
good  predictive  capabilities  of  the  similarity  models  for  correlations  between  the  modeled  and  exact  quantities.  They 
are  obtained  by  adding  an  eddy  viscosity  expression  to  a  similarity  model.  Mixed  models  were  found  to  provide  LES 
results  comparable  to  or  better  than  those  obtained  with  the  classical  Smagorinsky  model  with  the  additional  benefit 
of  accounting  for  the  backscatter  [27].  However,  because  the  results  do  not  dramatically  improve  and  the  dependence 
of  the  models  on  a  filter  introduces  an  additional  complication,  mixed  models  have  found  only  limited  acceptance. 

Due  to  the  diflBculties  encountered  by  the  above  procedures  in  describing  the  observed  properties  of  turbulence,  none 
of  them  has  been  accepted  as  a  fully  satisfactory  solution  to  the  problem  of  subgrid-scale  modeling.  This  motivates 
further  efforts  to  improve  the  status  of  SGS  modeling.  It  can  be  argued  that  improvements  in  modeling  can  only  be 
made  if  more  detailed  information  about  the  SGS  nonlinear  interactions  is  employed  than  information  available  to 
developers  of  the  older  models.  During  the  last  several  years  it  has  become  possible  to  investigate  in  great  detail  the 
nonlinear  interactions  in  turbulent  flows  at  low  Reynolds  numbers  using  direct  numerical  simulation  databases  [28-31] 
and  more  recently,  experimental  measurements  [24].  In  particular,  SGS  interactions  have  been  investigated  extensively 
for  isotropic  and  wall-bounded  turbulence  [14,16,32,33].  The  results  of  these  investigations  were  used  to  demonstrate 
that  many  of  the  observed  features  of  the  exact  SGS  interactions  can  be  inferred  from  the  dynamics  of  the  resolved 
scales  alone  [14,34]  (the  essential  goal  of  SGS  modeling)  using  procedures  entirely  different  from  the  traditional  SGS 
models.  These  procedures  are  mostly  of  theoretical  interest  since  they  have  not  been  implemented  in  actual,  time 
evolving  large  eddy  simulations.  However,  they  imply  a  possibility  of  improved  SGS  models  developed  by  using  the 
observed  properties  of  the  nonlinear  subgrid-scale  interactions.  The  purpose  of  this  paper  is  to  describe  and  evaluate 
a  new  approach  to  SGS  modeling  which  is  based  largely  on  concepts  derived  from  such  detailed  analyses  of  nonlinear 
interactions  in  a  variety  of  turbulent  flows. 


II.  VELOCITY  ESTIMATION  CONCEPT 

The  eddy  viscosity  models  and  the  similarity  models  illustrate  two  different  approaches  to  SGS  modeling.  The  eddy 
viscosity  models  assume  that  the  functional  form  of  the  SGS  stress  tensor  is  proportional  to  the  resolved  rate-of-strain 
and  attempt  to  provide  expressions  for  the  proportionality  coefficient,  the  eddy  viscosity.  In  the  similarity  approach, 
the  primary  modeled  quantity  is  the  full  velocity  field  and  no  functional  form  of  the  SGS  stress  is  assumed.  After  a 
model  for  the  full  velocity  field  is  introduced,  the  stress  tensor  is  computed  directly  from  the  definition  (5)  using  an 
appropriate  filter.  The  idea  of  estimating  the  unknown  velocity  field  from  the  resolved  one  for  the  purpose  of  SGS 
modeling  has  been  used  previously,  in  one  form  or  another,  by  several  authors.  For  example,  the  Bardina  model  is 
obtained  formally  by  approximating  the  full  velocity  by  the  filtered  velocity,  Ui  w  w,- ,  and  computing  according  to 
Eq.  (5).  Liu  et  al.  [24]  make  the  same  approximation  for  the  velocity,  but  use  a  wider  filter  to  compute  the  resolved 
stress  tensor  (9)  which  is  then  assumed  to  be  proportional  to  the  actual  SGS  stress  tensor  Shah  and  Ferziger  [35] 
have  recently  developed  a  procedure  to  compute  a  higher  order  approximation  to  the  full  velocity  field  than  the 
velocity  u,-  used  in  traditional  similarity  models.  The  approximated  velocity  and  a  modified  filter  definition  were  used 
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to  compute  the  SGS  stress  tensor  which  provided  good  results  in  LES  of  turbulent  channel  flow.  McDonough  proposed 
a  so-called  additive  turbulent  decomposition  [36]  for  simulations  of  turbulence  where  separate  equations  are  written  for 
large  and  small  scales  in  the  decomposition.  The  equations  for  large  (small)  scales  were  then  solved  with  a  proposed 
model  of  small  (large)  scales.  Jimenez  [37]  successfully  simulated  inertial  range  dynamics  using  a  numerical  scheme 
which  enforced  for  a  range  of  small  scales  the  4“®/^  energy  spectrum.  In  the  work  of  Kerr  et  al.  [14],  the  authors 
modeled  unresolved  subgrid-scale  vorticity  by  calculating  properly  normalized  vorticity  production  by  the  resolved 
scales  in  a  limited  range  of  wavenumbers  outside  the  resolved  range.  The  modeled  subgrid-scale  vorticity  was  used  in 
a  prior*  tests  to  compute  the  SGS  quantities  which  were  found  to  compare  very  favorably  with  the  exact  quantities. 
This  also  implies  that  the  nonlinear  interactions  among  resolved  scales  provide  significant  amounts  of  information 
about  the  largest  subgrid-scales.  This  conclusion  is  quite  consistent  with  the  classical  picture  of  turbulence,  where 
the  smaller  scales  are  generated  by  the  nonlinear  interactions  among  the  larger  scales. 

Despite  the  shortcomings  of  specific  methods  used  to  estimate  the  full  velocity  field,  the  above  results  taken  together 
suggest  that  this  approach  to  SGS  modeling  may  become  viable  if  better  estimation  procedures  are  developed.  In  this 
paper,  we  investigate  a  procedure  which  develops  a  significantly  better  model  for  the  velocity  field  than  the  above 
approximations.  The  procedure  explicitly  uses  properties  of  the  filtering  operation  and  properties  of  SGS  interactions 
established  in  previous  investigations. 

An  obvious  application  of  the  properties  of  a  filtering  operation  is  a  requirement  that  the  estimated  velocity  after 
filtering  with  a  given  filter  G  used  in  a  LES  is  equal  (or  close)  to  the  resolved  field  u,*  on  the  resolved  mesh.  This 
requirement  will  be  used  explicitly  in  the  next  section. 

The  properties  of  SGS  interactions  described  below  lead  to  additional  conditions  imposed  on  the  estimated  velocity. 
Subgrid-scale  energy  transfer  in  isotropic  turbulence  at  large  Reynolds  numbers  has  been  investigated  extensively  in 
the  framework  of  analytical  theories  of  turbulence  [5],  for  LES  databases  [31,32],  and  at  small  Reynolds  numbers 
using  well  resolved  DNS  fields  [16,32].  If  a  sharp  spectral  filter  with  a  cutoff  wavenumber,  kcj  is  used,  the  spectral 
eddy  viscosities  in  all  cases  exhibit  a  strong  cusp  in  the  vicinity  of  the  cutoff  wave  number.  The  cusp  behavior  reflects 
significant  SGS  energy  transfer  associated  with  modes  in  the  vicinity  of  the  cutoff.  Kraichnan’s  einalysis  [5]  for  an 
infinite  inertial  range  shows  that  close  to  75%  of  the  SGS  energy  transfer  is  from  the  range  of  the  resolved  scales 
0.5fcc  <  k  <  kc-  In  DNS  data  for  low  Reynolds  number  flows,  this  range  is  responsible  for  almost  the  entire  SGS  energy 
transfer  [16,32].  Additional  analysis  of  the  DNS  data  also  revealed  that  the  observed  SGS  energy  transfer  is  caused 
almost  exclusively  by  interactions  of  the  resolved  scales  with  a  limited  range  of  subgrid-scales  with  wavenumbers 
kc  <  k  <  2kc,  and  the  nonlocal,  eddy- viscosity  type  energy  transfer  between  resolved  modes  and  modes  k  >  2kc  was 
not  observed  [16].  At  high  Reynolds  numbers,  Kraichnan^s  analysis  implies  that  the  nonlocal  transfer  is  not  negligible 
and  may  contribute  about  25%  to  the  net  SGS  transfer.  These  results  strongly  suggest  that  the  SGS  energy  transfer  is 
dominated  by  energy  exchanges  among  resolved  and  unresolved  scales  from  the  vicinity  of  the  cutoff.  Similar  analyses 
have  been  performed  for  turbulent  wall-bounded  flows  (channel  flow  and  Rayleigh-Benard  convection)  and  the  same 
conclusions  have  been  reached  [33].  Also,  the  analysis  of  the  experimental  data  by  Liu  et  al.  [24]  identified  interactions 
between  neighboring  scales  as  dominant  in  the  SGS  dynamics.  A  major  conclusion  from  these  investigations  is  that 
the  nonlinear  dynamics  of  the  resolved  modes  with  wave  numbers  k  <  kc  are  governed  almost  exclusively  by  their 
interactions  with  a  limited  range  of  modes  with  wave  numbers  not  exceeding  2kc  and  that  much  smaller  scales  have  a 
negligible  effect  on  the  resolved  ones.  This  conclusion  is  entirely  consistent  with  the  classical  picture  of  the  turbulent 
cascade  [38],  where  the  large  scales  of  a  turbulent  flow  determine  the  energy  flux  down  the  spectrum  and  the  very 
small  scales  play  an  entirely  passive  role  by  adjusting  themselves  in  such  a  way  as  to  accommodate  this  energy  flux 
prescribed  by  the  large  scales.  In  the  context  of  SGS  modeling,  this  result  implies  that  in  the  definition  of  the  SGS 
stress  tensor  (5)  the  full  velocity  field  may  be  replaced  by  a  field  truncated  to  wavenumbers  k  <  2fcc,  Le.  to  scales 
not  smaller  than  about  half  of  the  smallest  resolved  scale.  More  generally,  consider  a  filtering  operation  (3)  with  two 
filter  widths  A  and  A/2  denoted  by  an  overbar  and  a  tilde,  respectively.  The  above  properties  of  the  SGS  dynamics 
imply  the  following  approximation 


Tij  —  UiUj  —  Ui  ILj  «  UiUj  —  Ui  Uj.  (10) 

Of  course,  the  approximation  in  its  present  form  does  not  provide  an  expression  for  in  terms  of  the  resolved 
quantities,  because  «,•  still  depends  on  the  unresolved  subgrid-scales  in  the  range  between  A/2  and  A.  However,  the 
approximation  (10)  simplifies  the  task  of  modeling  the  velocity  field,  because  for  the  purpose  of  SGS  modeling  its 
dependence  on  scales  smaller  than  A/2  can  be  neglected. 

Another  previously  mentioned  property  of  the  SGS  interactions  is  the  observation  that,  because  of  the  turbulent 
cascade,  the  subgrid-scales  in  the  range  between  A/2  and  A  must  be  strongly  influenced  by  the  resolved  scales  larger 
than  A.  In  particular,  during  the  initial  stages  of  the  evolution  of  a  turbulent  flow  the  smaller  scales  are  generated  by 
the  larger  eddies  and  the  phases  of  the  subgrid-scales  are  determined  by  the  nonlinear  interactions  among  the  large 
sc2Jes. 
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Using  the  qualitative  properties  of  filtering  and  SGS  interactions  described  above,  we  develop  in  the  subsequent 
sections  a  specific  method  for  estimating  the  velocity  «,*.  The  expectation  is  that  it  should  be  possible  to  develop  a 
good  model  for  the  velocity  «,•  by  accounting  for  the  physics  of  the  nonlinear  interactions  and  that  the  procedure  will 
result  in  a  better  model  for  the  full  SGS  stress  tensor  than  could  be  obtained  with  traditional  approaches  to  SGS 
modeling.  There  are  two  important  theoretical  advantages  of  such  an  approach.  First,  the  definition  (10)  guarantees 
that  the  model  based  on  the  estimated  velocity  is  Galilean  invariant  [39].  Second,  all  of  the  SGS  stress  components 
will  have  correct  near-wall  behavior,  because  the  estimated  velocity  5,*  is  approximated  to  the  lowest  order  by  the 
filtered  velocity  u,*  for  which  the  correct  wall  behavior  is  enforced  by  the  LES  equations. 


III.  A  VELOCITY  ESTIMATION  PROCEDURE 

Consider  a  velocity  field  Ui{x,t)  which  is  fully  resolved  in  a  particular  direction,  say  a:,  on  M  equally  spaced  grid 
points 


Xm  =  mS,  (m  =  0,l,...,Af-l)  (11) 

where  S  Lx/M  is  a.  mesh  size.  Such  a  field  could  be  obtained  in  adequately  resolved  direct  numerical  simulations. 
A  collocation  pseudo-spectral  method  with  Fourier  expansions  is  used  later  in  the  numerical  implementation  of  the 
model,  therefore  the  velocity  u,*  for  all  x  will  be  represented  by  a  trigonometric  polynomial  of  order  M  (Canuto  et 
al.  [40]) 


M/2 

^  Cmexp{ikmx), 
m=— Af/2 


(12) 


where  the  discrete  wavenumbers  are  km  =  27rm/Lx  =  mAk  and  the  dependence  of  on  other  variables  is  neglected 
in  the  notation.  The  superscript  M  signifies  the  M  mode  Fourier  representation  of  the  field  and  i  in  the  argument  of 
the  exponent  is  ^/—l.  The  Fourier  coefficients  Cm  in  (12)  are  expressed  in  terms  of  the  nodal  values  of  the  velocity 


M-l 


Cm  =  «.(*p)exp(-i*miCp), 


p=0 


(13) 


where  d  =  1  for  all  values  of  m  except  m  =  ±M  when  d  —  2.  The  velocity  field  «j'*^^(a:)  is  a  real  function  of  x,  thus 
the  complex  Fourier  coeflScients  must  satisfy  the  relation  =  c-m^  Equation  (13)  also  reveals  that  Cm  is  real  for 
m  =  0  and  m  =  ±iV.  All  of  the  coefficients  Cm  are  determined  by  M  real  numbers,  i.e.  the  same  number  of  nodal 
values  in  physical  space. 

In  Fourier  representation,  the  top  hat  filter  is  particularly  easy  to  implement.  The  filter  function  G  for  a  top  hat 
filter  with  a  width  A  is 


G(x, 


X 

A 


0 


if  \x  —  x'l  <  A/2 
otherwise. 


Integrating  (12)  with  the  filter  (14)  we  obtain 


M/2 

E 

m=— Af/2 


Cm 


sin(femA/2)' 

kmA/2  . 


exp(fAyy,x). 


(14) 


(15) 


The  filtered  field  is  smooth  over  the  length  scale  A  and  for  the  purpose  of  LES  it  can  be  accurately  represented  by 
sampling  it  on  a  courser  mesh  with  the  mesh  size  A.  It  is  convenient  to  assume  that  A  is  a  multiple  of  the  smaller 
mesh  size  A  =  rJ,  such  that  the  field  (15)  is  sampled  on  the  subset  of  2N  =  M/r  grid  points  from  the  fine  mesh. 
If  only  the  sampled  values  of  (15)  on  the  mesh  points  xj  =  jA,  {j  =  0,1,...,  2N  —  1)  are  used,  the  2N  mode  Fourier 
representation  of  U,*  is 


N 


u^^^\x)=  X,  aiexp{ikix), 


l=-N 


(16) 
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where  ki  =  2irlfLx  and 


-  2iV-l 

^  (17) 

i=0 

Since  (16)  is  the  same  quantity  as  the  velocity  Ui  in  actual  large  eddy  simulations,  the  superscript  2N  in  (16)  will 
often  be  omitted. 

The  maximum  wavenumber  in  the  representation  (16)  is  NAk.  According  to  the  discussion  in  the  previous  section, 
we  will  seek  a  model  velocity  u,-  which  contains  scales  by  at  most  a  factor  of  two  smaller  than  the  smallest  resolved 
scale.  In  Fourier  representation,  it  is  equivalent  to  doubling  the  maximum  wavenumber  which  requires  an  increase  to 
4N  modes  in  the  expansion 


2JV 

^  bnexp{iknx).  (18) 

n=-2JV 

The  modeling  procedure  attempts  to  express  4N  unknown  coefficients  in  (18)  in  terms  of  2N  known  coeflScients  a/ 
in  the  expansion  (16).  Such  a  procedure  will  not  be  unique,  but  if  the  observed  physics  of  the  nonlinear  interactions 
are  accounted  for,  good  models  for  a  SGS  stress  tensor  may  be  expected. 

We  propose  a  procedure  which  consists  of  two  steps,  one  which  is  purely  kinematical  and  the  other  which  accounts 
for  the  dynamics  of  SGS  nonlinear  interactions.  In  the  kinematic  step,  we  first  require  that  filtered  with  a  top 
hat  filter  with  the  width  A  is  equal  to  the  known  values  of  on  the  resolved  mesh 

U  =  0,1,...,2N  -  1).  (19) 

This  equation  reflects  the  compatibility  condition  between  the  resolved  velocity  and  the  filtered  estimated  velocity 
and  it  constitutes  2N  equations  relating  coefficients  6„  and  a/ .  The  additional  2N  equations  are  obtained  by  requiring 
that  and  filtered  with  a  wider  filter  (denoted  by  a  hat)  A  =  2A  are  equal  on  the  resolved  mesh 

{j  =  0,1,...,2N  -  1).  (20) 

The  use  of  the  wider  filter  in  the  above  condition  is  reminiscent  of  the  test  filtering  in  the  dynamic  model  [7].  While 
reasonable,  the  condition  (20)  is  less  justified  than  the  condition  (19)  since  it  does  not  reflect  the  straightforward 
properties  of  filtering  as  does  Eq.  (19).  Other  ways  of  obtaining  the  remaining  2N  equations  may  be  considered. 
For  example,  the  relation  (19)  may  be  enforced  on  points  interpolated  between  resolved  grid  points  or  the  coefficients 

bk  may  be  chosen  to  provide  the  optimal  approximation  to  «P^^(aT)  by  up^^(x)  (both  depending  on  the  continuous 
variable  x).  We  are  planning  to  explore  such  methods  in  future  work. 

Introducing  the  notation 


, _  sin(fcA/2) 
kA/2  ’ 

(21) 

_  sm(*A/2) 

^  kA/2  ’ 

(22) 

equations  (19)  and  (20)  are  rewritten  as  follows 


2N  .  N  ,  . 

^  a{kn)bn  exp(i^)  =  ^  a/  exp(i^),  j  =  0,1, . .  .,2N  - 

n--2N  l=-N  ^ 


(23) 


2N 


N 


X)  7(fcn)frn exp(i^)  =  ^  a,7(fe,)exp(j^),  J  =  0,l,...,2i\r-1 


N 

n=-2N 

Using  the  summation  formula  (see  Canuto  et  al,  [40]) 


N 


(24) 
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1  f  1  if p  =  2jVm,  m  =  0,±1,±2, ... 

^  =  {  0  otherwise. 

i=o  '• 

we  derive  a  set  of  equations  for  the  coefficients  b„.  Prom  equation  (23) 

ot(ki)bi  +  a(ki+2N)/>t+2N  =  oi,  —N  <l  <0 

a{ki)bi  +  Q{ki-.2N)bi-2N  =  ai,  0  <l  <  N, 

and  from  equation  (24) 

'yiki)bi  +  'y{ki+2N)bi+2N  =  'y{ki)ai,  —N  <  I  <  0, 


y(ki)bi  +  7(ki-2Jv)bi-2N  =  7(^/)<*»>  0  <  I  <  N. 

Solving  the  above  equations,  we  get  explicit  expressions  for  the  coefficients 

.  _  7(ki-2N)/7(ki)-a(ki-2N)  _  n  ^  ^  nr 

‘  a(ki)j(ki^2N)/7(ki) -0!(ki-2nr) 


(25) 

(26) 

(27) 

(28) 

(29) 

(30) 


“I 

In  special  cases  I  =  and  /  =  0,  we  obtain 

a(k. 


ine  loiiowmg  equations 


t(kr/)brr 


■  0‘{k-N)b-N  =  dN, 


7{kN)bN  +7(A:_w)6-Ar  =  7{kN)aN, 


(31) 

(32) 

(33) 


a(k2N)b2N  +  Oc{k^2N)b-2N  +  Ol(ko)bo  =  Uo, 


(34) 


7(k2N)b2N  +  7{k-2N)b-2N  +  7iko)bo  =  7iko)ao. 
Using  the  facts  that  a{ko)  =  7(^0)  =  1,  7(fc±2jv)  =  «(*±27v)  =  0,  and  7(Ar±jv)  =  0, 


(35) 

(36) 


and  the  coefficients  b2N  and  b^2N  can  be  set  to  zero.  It  can  be  shown  that  the  same  values  of  the  coefficients  6/ 
are  obtained  for  any  A  which  is  an  even  multiple  of  A.  Such  a  partial  independence  of  the  modeled  field  «,•  on  the 
secondary  filtering  operation  is  an  attractive  feature  of  the  modeling  procedure. 

The  above  method  provides  the  velocity  field  (18)  using  only  properties  of  the  filtering  operation  by  enforcing  the 
condition  that  the  filtered  modeled  field  is  equal  to  the  known  LES  field  (16)  on  a  selected  number  of  mesh  points. 
The  procedure  does  not  contain  any  information  about  the  fact  that  the  small  scales  in  such  a  field  are  dynamically 
coupled  to  large  scales  through  nonlinear  interactions.  In  order  to  account  for  such  couplings,  the  coefficients  bi  must 
be  modified.  According  to  the  discussion  in  section  II,  it  may  be  expected  that  the  smaller  scales  are  generated  by 
the  nonlinear  interactions  between  large  eddies.  In  spectral  methods,  the  large  scales  with  wavenumbers  k  <  kc  will 
generate  modes  with  wavenumbers  kc  <  k  <  2kc  through  the  nonlinear  term  in  the  Navier-Stokes  equations,  even  if 
initially  the  domain  kc  <  k  <  2k c  is  empty.  We  use  the  velocity  field  (18)  truncated  spectrally  to  k-N  <  k  <  kpfy  i.e. 


(*)  = 


N 

6„exp(»fc„x), 

n=-N 


(37) 


and  compute  for  this  field  the  nonlinear  term  for  wavenumbers  \k\  >  k^.  The  complex  coefficients  6^  for  |n|  >  iV  in 
the  original  expansion  (18)  are  modified  by  setting  their  phases  equal  to  the  phases  of  the  computed  nonlinear  term, 
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while  keeping  their  amplitudes  unchanged.  This  way,  the  highest  wavenumber  modes  in  the  Fourier  expansion  will 
carry  limited  information  about  the  process  of  generation  of  smaJl  scales  through  nonlinear  interactions.  The  actual 
phases  and  amplitudes  of  high  wavenumber  modes  are  generated  by  a  much  larger  number  of  interactions  acting  over 
at  least  one  large  eddy  turnover  time.  In  the  above  procedure,  this  process  is  simplified  by  assuming  that  only  a  single 
act  of  large  scale  interactions  determines  the  phases  of  subgrid-scales  and  their  amplitudes  can  be  determined  by  the 
filtering  compatibility  conditions  rather  than  by  time  evolution.  The  validity  of  these  assumptions  can  only  be  tested 
by  using  this  procedure  in  practice. 

In  summary,  the  procedure  described  above  allows  us  to  estimate  the  velocity  field  with  subgrid-scales  (18)  using 
only  information  about  the  filtered,  LES  field  (16).  The  only  free  parameter  is  the  additional  filter  width  A,  but  the 
estimated  velocity  is  the  same  if  A  is  any  even  multiple  of  the  primary  filter  width  A.  A  closed  form  expression  for 
the  estimated  velocity  cannot  be  easily  given,  because  of  the  modification  of  the  subgrid-scale  phases  by  the  nonlinear 
term.  The  generalization  of  the  estimation  procedure  to  more  dimensions  is  straightforward  if  Fourier  decomposition 
in  other  directions  can  be  used.  For  two  and  three  dimensional  fields,  the  kinematic  part  of  the  procedure  is  applied 
sequentially  in  each  direction,  followed  by  a  phase  correction  on  such  an  expanded  field.  Only  cases  allowing  Fourier 
expansions  will  be  considered  in  this  paper. 

As  already  noted,  the  procedure  is  not  unique  and  it  is  possible  that  it  could  be  improved  by  imposing  additional 
conditions  on  the  the  estimated  field  (18).  For  instance,  because  of  the  phase  correction  the  representation  (18)  no 
longer  satisfies  the  kinematic  conditions  (19)  and  (20)  exactly.  However,  these  conditions  are  satisfied  approximately 
because  the  large  scale  modes  and  the  amplitudes  of  the  small  scale  modes  are  unchanged.  The  current  procedure 
also  does  not  enforce  incompressibility  for  the  estimated  field  (18).  On  the  other  hand  it  must  be  remembered  that 
estimating  velocity  is  not  a  goal  in  itself,  but  an  intermediate  step  in  calculating  the  SGS  stress  tensor.  If  the 
estimation  procedure  properly  captures  the  kinematic  and  dynamic  relations  between  the  resolved  and  subgrid-scales, 
it  should  produce  a  good  approximation  for  the  SGS  stress  tensor. 


IV.  NUMERICAL  IMPLEMENTATION 

We  will  consider  turbulent  channel  flow  as  a  test  case,  because  of  the  availability  of  experimental  and  DNS  data 
for  comparisons.  The  fluid  is  contained  between  two  parallel  rigid  walls  separated  by  a  distance  Lz  in  the  vertical 
direction  z.  The  horizontal  dimensions  of  the  computational  domain  are  Lx  in  the  streamwise  direction  x  and  Ly 
in  the  spanwise  direction  y.  In  the  simulations,  no-slip  boundary  conditions  are  imposed  on  the  rigid  walls  and 
periodic  boundary  conditions  in  x  and  y.  We  employ  a  numerical  Navier-Stokes  code  developed  by  Chan  [41]  which 
uses  a  pseudo-spectral  numerical  method  for  spatial  discretization  [40]  with  Fourier  expansions  in  the  streamwise  and 
spanwise  directions,  and  Legendre  polynomials  in  the  wall  normal  direction.  The  equations  are  integrated  in  time 
using  a  second  order,  predictor-corrector  method  introduced  by  Kim  and  Moin  [42].  The  LES  equations  (1)  and 
(2)  are  nondimensionalized  by  the  channel  half  width  h  =  Lzl2  and  the  friction  velocity  «r  =  where  tq 

is  the  wall  shear  stress.  The  pressure  term  is  decomposed  into  the  fluctuating  component  and  the  mean  pressure 
gradient  driving  the  flow  which,  for  the  nondimensionalization  applied,  is  equal  to  —1  for  the  streamwise  direction. 
The  nondimensionalized  equations  are 


a  _  _ 


Rcr  dxjdxj  dx^ 


-Ti 


tji 


(38) 


where  Rcr  =  Urhjv  is  the  Reynolds  number.  The  filtering  in  (38)  and  (39)  is  applied  only  in  the  horizontal  directions 
with  the  top  hat  filter  with  widths  Ai  and  A2  in  x  and  j/,  respectively.  The  widths  are  chosen  to  be  equal  to 
the  respective  mesh  sizes  in  LES  for  those  directions.  Since  Fourier  expansions  are  implemented  in  x  and  j/,  a  two- 
dimensional  version  of  the  estimation  method  described  in  the  previous  section  is  used  to  model  the  SGS  stress  tensor, 
Tiy  For  comparison,  we  also  use  the  Smagorinsky  model  which  provides  the  following  expression  for  the  eddy  viscosity 
in  Eq.  (6) 


i/t  =  2(CsA)2|5|  (40) 

where  |5|  =  {2Sij  and  the  length  scale  A  is  taken  as  the  geometric  average  of  mesh  spacings  in  the  Cartesian 

directions  A  =  (AxAyAz)^/^.  The  value  of  the  Smagorinsky  constant  was  estimated  theoretically  by  Lilly  [43]  as 
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Cs  =  0.18,  but  we  use  Cs  =  0.1  which  is  the  preferred  value  for  wall-bounded  shear  flows  [25,44].  In  order  to 
control  excessive  dissipation  by  the  model  in  the  vicinity  of  walls,  the  length  scale  A  is  decreased  using  Van  Driest 
damping  [45]  as  proposed  by  Moin  and  Kim  [46]. 


A.  A  priori  analysis  of  DNS  data 

In  order  to  assess  the  validity  of  the  assumptions  made  in  developing  the  velocity  estimation  model,  the  modeled 
SGS  quantities  are  compared  in  this  section  with  the  exact  SGS  quantities  computed  for  well  resolved  DNS  databases 
obtained  in  pseudo-spectral,  Fourier-Chebyshev  simulations  of  turbulent  channel  flow  described  by  Gilbert  [47]  and 
Gilbert  and  Kleiser  [48].  The  flow  was  simulated  with  a  resolution  of  160^(horizontal)x  129  (vertical)  modes  in  a 
box  with  dimensions  =  2  (the  unit  of  length  is  the  channel  half  width,  h  =  Lzl2)^  Lx  =  S.Gtt,  and  Ly  —  I.Ott. 
Simulations  were  performed  enforcing  constant  mean  mass  flux  at  Reynolds  number  Re  =  5000  based  on  the  maximum 
velocity  of  a  laminar  flow  and  the  channel  half  width  ft.  The  corresponding  Reynolds  number  based  on  friction  velocity 
and  h  is  approximately  Rtr  =  210.  In  the  spanwise  direction  a  symmetry  condition  at  y  =  Ly/2  was  imposed, 
motivated  by  the  symmetry  properties  of  the  transitional  flow  preceding  the  fully  turbulent  stage  of  the  simulation.  A 
comparison  with  the  DNS  results  of  Kim  et  al,  [49]  and  experimental  results  of  Nishino  and  Kasagi  [50]  revealed  that 
the  influence  of  the  symmetry  assumption  on  horizontally  averaged  statistical  properties  of  turbulence  was  negligible. 
The  simulations  were  run  until  turbulence  reached  a  statistically  steady  state,  and  after  that  they  were  continued 
for  long  enough  times  to  collect  turbulence  statistics.  For  a  priori  analyses  only  instantaneous  fields  were  used  with 
averaging  over  planes.  Averaging  procedures  will  be  denoted  by  brackets  <  ...  >  to  distinguish  them  from  the  filtering 
operation  which  is  denoted  by  the  overbar. 

The  artificial  LES  velocity  field  was  obtained  from  the  DNS  velocity  as  follows.  First,  in  the  wall  normal  direction 
the  velocity  was  interpolated  from  129  Chebyshev  grid  points  to  65  Legendre-Gauss-Lobatto  points  on  the  interval 
— 1  <  2:  <  +1.  Second,  the  top  hat  filtering  in  the  form  of  (15)  was  applied  with  Ai  =  Lar/32  and  A2  =  Ly/Z2  in  the 
horizontal  directions  x  and  y,  respectively.  The  continuous  x  and  y  representation  (15)  was  subsequently  sampled  on 
32  X  32  grid  points  uniformly  spaced  in  each  direction.  The  resulting  LES  field  is  thus  represented  on  a  mesh  with 
32  X  32  X  65  points.  The  LES  velocity  was  used  to  compute  SGS  quantities  using  both  the  Smagorinsky  model  and  the 
velocity  modeled  by  the  estimation  procedure  described  in  section  III.  The  exact  SGS  quantities  were  calculated  from 
definitions  using  the  full  DNS  velocity  and  top  hat  filtering  in  x  and  y  with  the  respective  filter  widths  Ai  =  La:/32 
and  A2  =  Ly/32. 

In  addition  to  the  SGS  stress  tensor  several  other  SGS  quantities  were  computed:  the  subgrid-scale  nonlinear 
term  (or  SGS  force) 


(41) 

the  SGS  energy  transfer 

(42) 

and  the  SGS  dissipation 

esGs(x)  =  Tij{x)Sij{x). 

(43) 

The  primary  quantity  that  must  be  properly  predicted  by  a  SGS  model  is  the  SGS  stress  tensor.  All  of  the  SGS 
stress  tensor  components  were  calculated  exactly  from  the  definition  Eq.  (5),  and  using  the  Smagorinsky  and  the 
velocity  estimation  models.  In  figure  1,  we  plot  results  for  four  nonvanishing  components  of  the  tensor  averaged  over 
the  horizontal  planes.  The  diagonal  components  are  modified  by  subtracting  one  third  of  the  trace  of  the  tensor, 
because  this  is  the  quantity  modeled  by  the  Smagorinsky  expression  (6).  The  Smagorinsky  model  predicts  essentially 
a  zero  value  for  the  plane  averaged  diagonal  components  of  the  stress.  The  modeled  SGS  shear  stress  component  ri3 
is  non-zero,  but  it  is  significantly  under  predicted  throughout  the  channel.  Note,  however,  that  the  shear  stress  slope 
is  predicted  very  well  at  the  walls.  For  the  diagonal  components,  the  estimation  model  properly  predicts  trends  with  z 
and  locations  of  peaks,  though  the  peak  values  are  underestimated  by  up  to  50%.  The  estimation  model  also  provides 
qualitatively  better  results  for  ri3  across  the  channel,  but  the  slopes  at  the  walls  are  slightly  under  predicted. 

The  plane  averages  of  the  SGS  force  components  (41)  are  plotted  in  figure  2.  Again,  for  comparison  with  the 
Smagorinsky  model  the  expression  (41)  is  computed  with  the  diagonal  stress  components  modified  as  above.  The 
Smagorinsky  model  predicts  nonzero  values  only  for  the  streamwise  force  component  ,  Due  to  the  previously 

noted  favorable  prediction  for  the  wall  slope  of  ri3,  the  SGS  force  is  in  very  good  agreement  with  the  exact  quantity  in 
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the  immediate  vicinity  of  the  wall,  but  the  negative  peak  in  Fig.  2(a)  is  too  narrow.  The  estimation  model  captures 
correctly  the  locations  and  widths  of  the  peaks  in  the  exact  SGS  forces,  but  underestimates  their  extrema.  A  relatively 
good  prediction  of  the  estimation  procedure  for  the  wall  normal  force  implies  that  the  wall  normeJ  pressure 

gradient  may  be  equally  well  predicted  since  the  pressure  term  does  not  have  to  compensate  for  the  errors  in  the  wall 
normal  SGS  force.  If  the  Smagorinsky  model  is  used  the  pressure  in  Eq.  (1)  is  replaced  by  the  sum,  p+  l/S'Hck^ij, 
and  the  pressure  itself  cannot  be  computed,  because  the  model  does  not  provide  an  expression  for 

The  plane  averaged  SGS  dissipation  is  plotted  in  figure  3(a).  Both  models  give  correct  qualitative  trends  and  the 
Smagorinsky  model  prediction  is  very  accurate  in  the  wall  region  where  the  estimation  model  under  predicts  the  excict 
quantity.  The  favorable  performance  of  the  Smagorinsky  model  is  not  surprising  since  it  is  designed  to  account  for 
the  dissipative  features  of  turbulence  and,  moreover,  we  use  the  value  of  the  Smagorinsky  constant  optimized  for  the 
modeled  case.  In  figure  3(b),  we  plot  the  SGS  dissipation  decomposed  in  each  plane  into  forward  transfer  (negative 
values,  <  >)  and  inverse  transfer  or  backscatter  (positive  values,  <  €gQg  >).  The  actual  backscatter  is  a 

function  of  the  filter  and  is  fairly  small  for  the  top  hat  filter  used  here.  Nevertheless,  the  estimation  model  predicts 
backscatter  in  excellent  agreement  with  the  exact  result  while  the  Smagorinsky  model  is  purely  dissipative. 

We  have  also  calculated  correlation  coefficients  between  the  modeled  and  exact  fields  for  a  number  of  subgrid-scale 
quantities.  The  correlation  coefficient  between  two  scalar  fields  Ti(x)  and  T2(x)  is 


C(Ti,T2) 


<  ri(x)r2(x)  >  ~  <  ri(x)  ><  r2(x)  > 

^((Ti(x))2>  -  (Ti(x))y<(T2(x))2)  - 


(44) 


Here,  <  ...  >  denotes  averages  taken  over  all  mesh  points  (a  global  correlation  coefficient)  or  over  points  in  the 
horizontal  planes  (a  plane  correlation  coefficient).  In  the  latter  case  the  correlation  coefficient  is  a  function  of  the 
vertical  coordinate  z.  The  global  correlation  coefficients  are  summarized  in  Table  I. 

The  correlations  for  individual  components  of  the  SGS  stress  tensor  are  around  0.5  for  the  estimation  model,  with 
the  exception  of  the  T23  component,  and  are  much  larger  than  the  correlations  for  the  Smagorinsky  model.  Both 
models  give  similar  correlation  coefficients  for  the  streamwise  component  of  the  SGS  force  as  well  as  for  the  transfer 
and  dissipation.  However,  the  wall  normal  SGS  force  in  the  estimation  model  shows  a  much  higher  correlation  than 
the  force  in  the  Smagorinsky  model.  The  plane  correlation  coefficients  for  the  important  ria  component  is  shown  in 
figure  4.  In  the  midsection  of  the  channel  the  estimation  model  gives  very  high  values  of  the  correlation  coefficient, 
about  0.6,  comparable  to  values  observed  in  consistent  a  priori  tests  for  similarity  models  [24].  In  the  same  region, 
the  Smagorinsky  model  gives  much  lower  correlations,  about  0.2.  The  correlations  decrease  in  the  vicinity  of  the  walls 
and  even  become  negative  for  the  Smagorinsky  model.  The  low  correlations  in  the  vicinity  of  the  walls  may  suggest 
that  the  model  predictions  are  deficient  in  this  important  region  where  the  mean  stresses  reach  meiximum  values.  On 
the  other  hand,  it  is  interesting  to  no^  in  figure  5  that  the  estimation  model  predicts  the  cross-correlation  between 
ri3  and  the  rate  of  strain  component  S13  in  very  good  agreement  with  the  correlations  for  the  exact  quantities  even 
in  the  wall  region.  One  may  speculate  that  in  a  time  evolving  LES  any  inaccuracies  in  the  prediction  of  the  SGS 
stress  in  the  wall  region  will  be  compensated  by  the  fact  that  the  modeled  stress  is  correlated  with  other  physical 
quantities  in  a  manner  similar  to  the  actual  stress.  This  hypothesis  can  only  be  confirmed  by  performing  aotual  large 
eddy  simulations. 

In  order  to  assess  the  dependence  of  the  results  on  the  LES  resolution,  additional  a  priori  tests  were  also  performed 
using  DNS  fields  filtered  to  48  x  64  x  65  grid  points.  Improved  comparisons  with  the  exact  results  and  increased 
values  of  the  correlation  coefficients  were  observed  for  all  SGS  quantities.  These  results  are  not  reported  here  in 
detail,  because  the  actual  LES  results  presented  in  the  next  section  were  obtained  from  simulations  implementing  the 
lower  resolution,  32  x  32  x  65  points. 

The  main  conclusions  from  a  priori  tests  are  that  the  SGS  modeling  procedure  based  on  the  velocity  estimation 
predicts  most  of  the  SGS  quantities  in  good  agreement  with  the  exact  SGS  quantities  and  performs  much  better 
overall  than  the  classical  Smagorinsky  model  optimized  for  the  evaluated  flow.  In  order  to  fully  assess  its  utility, 
it  must  be  implemented  and  tested  in  time  evolving  large  eddy  simulations.  The  results  from  such  simulations  are 
presented  in  the  next  subsection. 


B-  Large  eddy  simulations 

The  Smagorinsky  and  the  velocity  estimation  model  were  implemented  in  a  computer  code  to  solve  Eqs.  (38) 
and  (39).  The  parameters  used  in  these  simulations  are  presented  in  Table  11.  The  computations  A,  SI,  and  El 
for  Rcr  =  180  are  standard  test  cases  for  LES,  because  of  the  availability  of  several  independent  numerical  and 
experimental  databases  for  turbulent  channel  flow  [47,49,50].  Case  A  is  an  under-resolved  direct  numerical  simulation 
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and  is  used  only  for  comparison  with  cases  SI  and  El  to  provide  a  qualitative  estimate  of  the  influence  of  a  SGS 
model  on  the  flow  evolution.  Cases  S2  and  E2  correspond  to  one  of  the  simulations  of  Piomelli  [51]  performed  with 
the  dynamic  model  which  compared  well  with  the  experimental  results  of  Wei  and  Willmarth  [52].  In  all  cases  the 
dimensions  of  the  computational  domain  were  equal  to  the  dimensions  used  in  the  DNS  of  channel  flow  by  Kim  et 
al.  [49]  {Lx  =  47r,  Ly  =  27r,  and  =  2),  and  comparable  to  those  used  by  Gilbert  [47].  For  the  higher  Reynolds 
number  cases  S2  and  E2,  the  horizontal  dimensions  are  larger  than  those  used  by  Piomelli  [51]  {Lx  =  2.57r,  Ly  =  0.57r). 
The  number  of  grid  points  is  typical  of  other  large  eddy  simulations  at  these  Reynolds  numbers.  In  cases  S2  and  E2, 
the  number  of  grid  points  is  the  same  as  in  Case  3  of  Piomelli  [51],  but  the  grid  resolution  is  lower,  because  of  the 
larger  size  of  the  computational  domain. 

All  of  the  cases  were  started  from  the  same  initial  condition  which  was  generated  in  the  following  manner.  The 
code  with  the  Smagorinsky  model  was  initialized  with  one  of  the  DNS  fields  of  Gilbert  [47]  and  run  until  a  statistically 
steady  state  was  reached.  The  restart  file  from  the  end  of  this  run  was  used  as  the  initial  condition  for  all  cases  which 
were  run  for  the  times  (nondimensionalized  by  Ur/h)  indicated  in  Table  II.  A  constant  time  step  At  =  0.001  was 
used  in  all  cases.  Cases  with  the  estimation  model  required  about  40%  more  time  per  time  step  than  cases  with  the 
Smagorinsky  model.  An  important  quantity  monitored  for  each  case  was  the  shear  velocity,  Ur ,  which  was  computed 
from  the  mean  velocity  in  the  simulations.  For  the  nondimensionalization  used,  Ur  must  be  equal  to  unity  at  steady 
state  and  in  each  case  Ur  was  always  within  5%  of  the  steady  state  value.  This  is  the  typical  accuracy  observed  in 
other  large  eddy  simulations  of  a  channel  flow  [53].  We  found  that  restarting  the  higher  Reynolds  number  case  E2  from 
the  lower  Reynolds  number  case  led  to  a  very  slow  evolution  of  the  mean  flow,  and  the  computation  required  more 
than  22  time  units  for  the  shear  velocity  to  approach  the  steady  state  value  within  5%.  The  run  time  was  shortened 
in  case  S2  by  setting  the  mean  velocity  in  the  initial  condition  to  the  value  estimated  for  the  final  steady  state.  The 
turbulence  statistics  were  collected  over  last  two  time  units  and  they  typically  involved  ten  independent  realizations. 
For  consistent  comparison  with  previous  results,  the  raw  data  from  the  current  simulations  were  renormalized  with 
Ur  computed  from  the  mean  velocity. 

All  of  the  SGS  quantities  computed  and  analyzed  a  priori  in  the  previous  section  were  also  calculated  in  the  LES 
cases  SI  and  El.  The  general  behavior  observed  in  the  LES  is  very  similar  to  that  observed  in  a  priori  analysis.  In 
figure  6,  we  plot  the  non-vanishing  components  of  the  SGS  stress  tensor  with  one  third  of  the  trace  of  the  tensor 
subtracted  for  the  diagonal  components.  As  in  the  a  priori  analysis,  the  diagonal  components  are  zero  for  the 
Smagorinsky  model.  All  of  the  stress  components  for  the  estimation  model  have  peak  values  larger  than  in  a  priori 
analysis  and  are  in  much  better  agreement  with  the  exact  quantities  than  could  be  expected  from  a  priori  results. 
Similar  behavior  is  also  observed  for  the  SGS  force  components  plotted  in  figure  (7).  The  results  obtained  in  case  SI 
are  consistent  with  a  priori  analysis  while  the  SGS  force  from  case  El  is  in  much  closer  agreement  with  the  exact 
quantity  for  both  force  components  plotted.  These  results  suggest  that  the  velocity  estimation  procedure  is  more 
robust  than  a  priori  analyses  indicate.  We  speculate  that  the  modeled  velocity  must  adjust  itself  to  more  appropriate 
values  through  repeated  application  of  the  estimation  procedure  over  many  time  steps. 

The  plane  averaged  SGS  dissipation  for  cases  SI  and  El  is  plotted  in  figure  8(a).  The  dissipation  peaks  are  higher 
in  case  SI  than  in  case  El,  but  the  difference  is  much  less  than  what  was  observed  in  a  priori  analysis  and  the  SGS 
dissipation  in  both  cases  agrees  quite  well  with  the  exact  quantity.  This  feature  of  the  model  is  quite  important, 
because  it  distinguishes  it  from  other  non-eddy  viscosity  models  which  are  known  to  seriously  under  predict  SGS 
dissipation.  In  figure  8(b),  we  plot  the  SGS  dissipation  decomposed  in  each  plane  into  forward  transfer  (negative 
values,  <  backscatter  (positive  values,  <  c^qs  >)•  The  forward  transfer  in  case  El  is  very  close  to  the 

exact  forward  transfer  and  to  the  total  SGS  dissipation  of  the  Smagorinsky  model.  However,  contrary  to  the  results 
of  a  priori  analysis  the  modeled  backscatter  is  now  larger  than  the  actual  backscatter.  The  values  of  the  forward 
transfer  and  the  backscatter  integrated  across  the  channel  for  all  of  the  LES  cases  and  the  exact  results  for  the  DNS 
data  are  shown  in  Table  III.  The  SGS  dissipation  results  for  case  El  are  important,  because  they  demonstrate  that 
LES  with  the  estimation  model  can  be  run  for  thousands  of  time  steps  in  a  numerically  stable  manner  with  significant 
local  backscatter.  This  distinguishes  the  estimation  model  from  the  dynamic  model  where  the  backscatter  inherent 
in  the  model  must  be  limited  by  additional  procedures  in  order  to  prevei^  numerical  instabilities. 

The  plane  correlations  between  na  and  the  rate  of  strain  component  5i3  in  LES  are  shown  in  figure  9.  A  slight 
improvement  in  the  prediction  of  the  correlation  peaks  for  the  estimation  model  is  observed.  The  good  quality  of  the 
comparison  with  the  correlation  coefficients  for  the  exact  quantities  implies  that  the  correct  correlations  between  t\z 
and  5i3  are  maintained  in  time  evolving  large  eddy  simulations  using  the  estimation  procedure. 

The  SGS  stresses  are  intermediate  quantities  modeled  to  predict  mean  and  fluctuating  turbulent  velocities  which 
are  the  primary  quantities  of  interest  in  actual  LES.  To  obtain  these  quantities  from  the  computational  data  of  cases 
SI,  El,  and  A,  averages  were  taken  over  ten  realizations  and  over  each  horizontal  plane.  The  exact  DNS  results 
were  also  horizontally  averaged,  but  only  four  statistically  independent  realizations  were  available.  In  figure  10,  we 
plot  mean  velocity  in  wall  coordinates  for  the  well  resolved  DNS  data  and  for  the  three  cases  A,  SI,  and  El.  The 
standard  linear  law  of  the  wall  and  the  logarithmic  law  with  the  additive  constant  5.5  are  also  included  in  the  plot. 
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In  the  low  resolution  DNS  case  A,  the  mean  velocity  is  under  predicted  in  the  log  region.  The  mean  velocity  in 
case  SI  slightly  underestimates  the  exact  DNS  result  up  to  =  30,  but  then  overestimates  it  in  the  log  region 
>  30.  In  case  El,  the  mean  velocity  follows  the  exact  DNS  result  very  closely  up  to  =  10  —  20,  and  then  slightly 
underestimates  it  throughout  the  rest  of  the  channel  >  20.  Overall,  in  all  of  the  cases  the  mean  velocity  profiles 
are  in  reasonable  agreement  with  the  exact  DNS  result  throughout  the  wall  region,  but  then  vary  in  the  log  region. 
In  particular,  case  El  reveals  improvement  over  the  under  resolved  DNS  case  A  and  lies  closer  to  the  exact  DNS 
results  than  the  Smagorinsky  case  S2.  The  plane  averaged  root-mean-square  (rms)  turbulent  velocities  are  plotted 
in  figure  11.  The  turbulent  velocity  in  LES  is  defined  as  it,"  =  u,-  —  [/i,  where  i/,-  is  the  plane  and  time  averaged 
mean  velocity  Ui  =<  u,*  >.  This  quantity  differs  from  the  standard  fluctuating  velocity  u/  =  u,—  <  u,-  >,  which 
includes  contributions  from  the  unknown  subgrid-scales.  Frequently,  in  comparing  LES  results  with  the  experimental 
data,  this  distinction  is  not  made,  because  experiments  only  provide  the  fluctuating  velocity  li,'.  However,  consistent 
comparisons  can  be  made  if  well  resolved  DNS  data  are  used.  In  Fig.  11,  we  compare  rms  velocities  predicted  by  the 
models  with  the  exact  DNS  results  calculated  for  ti,'  and  First  we  concentrate  on  the  consistent  comparison  with 
the  filtered  DNS  results  In  case  A,  we  observe  that  without  a  SGS  model  Vrms  and  Wrms  are  significantly  over 
predicted  compared  with  the  filtered  DNS  results  and  the  predictions  of  the  models.  This  points  to  an  appreciable 
effect  of  the  models  on  the  simulations.  However,  for  Urms  case  A  gives  results  intermediate  between  cases  SI  and  El. 
In  cases  SI  and  El,  the  peak  values  in  Urms  are  over  predicted  by  about  25%  and  the  Smagorinsky  model  also  over 
predicts  Urms  inside  channel.  The  estimation  model  gives  a  good  prediction  for  z*^  >  30.  For  the  other  two  velocity 
components,  both  models  over  predict  the  exact  rms  values  with  the  Smagorinsky  model  giving  slightly  better  results. 
As  expected,  the  rms  values  for  the  unfiltered  velocities  u,'  are  greater  than  those  of  the  filtered  velocities  ti,".  These 
results  also  illustrate  that  different  conclusions  about  the  performance  of  the  models  may  be  reached  if  only  rms 
values  of  the  unfiltered  velocities  ii,'  are  available.  Indeed,  in  such  a  case  the  Smagorinsky  model  gives  the  best  results 
for  Urms  and  the  new  model  is  better  for  the  other  two  velocity  components.  The  differences  between  filtered  and 
unfiltered  rms  velocities  will  be  even  larger  for  higher  Reynolds  numbers,  because  the  subgrid-scales  contain  a  larger 
fraction  of  the  total  turbulent  kinetic  energy.  Therefore,  comparisons  between  rms  velocities  in  LES  and  experiments 
must  be  interpreted  with  caution. 

The  main  conclusion  from  the  results  for  the  mean  and  fluctuating  turbulent  velocities  is  that  the  estimation  model, 
without  any  adjustable  constants,  gives  results  of  comparable  quality  to  the  results  obtained  with  the  Smagorinsky 
model  optimized  for  channel  flow  by  a  proper  choice  of  the  Smagorinsky  constant  and  the  wall  damping  function. 

The  effects  of  grid  resolution  were  investigated  by  changing  the  resolution  to  48  x  64  and  16  x  16  mesh  points 
in  the  horizontal  directions.  For  both  models,  the  increased  resolution  improved  predictions  of  the  rms  velocities, 
but  the  mean  velocity  was  slightly  over  predicted.  Nevertheless,  the  conclusions  previously  drawn  about  the  models 
remain  unchanged.  For  the  lower  resolution  of  16  x  16  points,  the  predictions  deteriorated  markedly,  especially  for 
the  streamwise  rms  velocity.  Inspection  of  the  horizontal  energy  spectra  in  the  fully  resolved  DNS  revealed  that  this 
resolution  corresponded  to  a  truncation  of  spectra  just  before  the  energy  peaks  and  that  the  resolution  of  32  x  32 
points  was  adequate  to  capture  the  peaks.  These  observations  suggest  that  the  minimum  resolution  requirements  for 
reasonable  LES  results  is  that  the  energy  containing  range  is  resolved  in  the  simulations. 

Cases  S2  and  E2  were  done  to  assess  the  performance  of  the  estimation  procedure  for  high  Reynolds  numbers  flows. 
The  results  of  the  large  eddy  simulations  of  Piomelli  [51]  performed  with  the  dynamic  model  at  the  same  Reynolds 
number  serve  as  a  useful  benchmark  dataset.  However,  there  are  several  differences  between  his  simulations  and  cases 
S2  and  E2  which  need  to  be  pointed  out.  Piomelli  used  the  spectral  cutoff  filter  while  we  used  the  tophat  filter,  and 
his  simulations  were  better  resolved  with  the  same  number  of  grid  points,  because  his  computational  domain  was 
smaller.  Also,  the  Chebyshev  mesh  used  by  Piomelli  for  the  normal  direction  had  better  resolution  in  the  wall  region 
than  the  Legendre  mesh  employed  in  our  simulations. 

In  figure  12  we  present  mean  velocity  profiles  in  wall  units  for  cases  S2  and  E2  and  for  Piomelli ’s  [51]  case  3. 
In  case  S2  with  the  Smagorinsky  model,  the  mean  velocity  approaches  the  logarithmic  curve  away  from  the  walls, 
but  it  is  significantly  underestimated  for  10  <  z"*"  <  100.  It  appears  that  this  problem  is  related  to  the  particular 
value  Cs  =  0.1  of  the  Smagorinsky  constant  chosen  here.  Piomelli  estimates  the  effective  Smagorinsky  constant  for 
his  simulations  as  0.06.  Shah  and  Ferziger  [35]  used  Cs  =  0.065  which  resulted  in  the  mean  flow  overshooting  the 
logarithmic  curve  and  Piomelli’s  data  for  z"*"  >  10.  The  mean  velocity  in  case  E2  departs  from  the  Piomelli’s  result 
around  z"*"  «  20  and  lies  above  the  logarithmic  curve  throughout  most  of  the  channel  width. 

The  rms  turbulent  velocities  are  plotted  in  figure  13.  The  curves  for  the  Smagorinsky  model  do  not  exhibit  gradients 
at  the  wall  as  steep  as  the  other  two  curves  and  they  generally  lie  below  Piomelli’s  data.  The  exception  is  in  Urms 
where  the  peak  in  the  vicinity  of  the  wall  is  broader  for  case  S2  than  for  the  other  two  cases,  and  the  Smagorinsky 
curve  lies  above  Piomelli’s  results  for  —0.96  <  z  <  —0.7.  For  a  lower  value  of  the  Smagorinsky  constant,  Cs  —  0.065, 
Shah  and  Ferziger  [35]  obtained  the  location  of  the  peak  in  better  agreement  with  Piomelli’s  results,  but  the  peak 
value  was  over  predicted  (urma/^T  ^  3.5),  and  for  \z\  <  0.9,  the  Urms  results  were  consistently  below  Piomelli’s  curve. 
In  case  E2,  the  rms  velocities  outside  of  the  wall  boundary  layers  are  consistently  below  the  values  predicted  by  the 
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Smagorinsky  and  the  dynamic  model.  As  noted  previously,  the  performance  of  the  models  is  diflScult  to  judge  based 
on  the  rms  velocities  because  of  differences  between  the  respective  LES  and  actual  turbulent  velocities  ti,"  and  u/.  In 
fcict,  the  rms  velocity  results  for  the  low  Reynolds  number  cases  suggest  that  good  SGS  models  should  predict  lower 
rms  velocities  in  LES  than  measured  experimentally.  This  is  the  case  for  run  E2,  because  Piomelli’s  [51]  data  match 
or  lie  below  the  experimental  results  of  Wei  and  Willmarth  [52].  However,  without  more  detailed  information  about 
the  scale  decomposition  of  the  energy  in  experiments  it  is  not  possible  to  decide  which  model  gives  the  best  results. 
In  the  immediate  vicinity  of  the  wall,  the  curves  from  case  E2  exhibit  the  same  steep  gradients  observed  in  Piomelli’s 
data  and  similar  locations  of  the  peak  values,  but  the  peak  values  themselves  are  not  the  same.  In  particular,  the  peak 
values  of  Urma  in  case  E2  are  larger  than  in  the  other  two  cases.  This  behavior  seems  to  suggest  that  the  estimation 
procedure  provides  less  dissipation  in  the  wall  region  and  more  dissipation  in  the  mid-channel  than  the  Smagorinsky 
and  the  dynamic  models.  This  conclusion  is  only  partially  supported  by  the  SGS  dissipation  results  plotted  in  figure 
14.  It  is  seen  that  in  the  wall  region  —1.0  <  z  <  —0.9  the  Smagorinsky  model  indeed  predicts  substantially  higher 
dissipation  than  the  estimation  model.  Outside  of  the  wall  region  the  forward  transfer  predicted  by  the  estimation 
model  is  the  same  as  that  predicted  by  the  Smagorinsky  model,  but  the  total  dissipation  is  lower,  because  of  the 
non-vanishing  backscatter  in  the  former  case.  The  lower  SGS  dissipation  in  case  E2  may  thus  explain  higher  Urms 
peaks  than  observed  in  case  S2,  but  lower  values  of  all  rms  velocities  in  the  mid-channel  are  not  consistent  with  such 
a  simple  energetic  argument. 

The  decomposition  of  the  SGS  dissipation  into  positive  and  negative  components  in  Fig.  14  reveals  substantial 
values  of  backscatter  in  case  E2.  According  to  Table  III,  the  backscatter  equals  60%  of  the  net  SGS  dissipation.  The 
local  physical  space  backscatter  appears  naturally  in  the  procedure  and  did  not  lead  to  any  numerical  instabilities  in 
simulations  extending  for  thousands  of  time  steps.  The  net  SGS  dissipation  in  case  S2  is  much  larger  than  in  case 
E2.  Decreasing  the  Smagorinsky  constant  from  0.1  to  0.06  in  the  definition  of  the  SGS  dissipation  would  reduce  the 
dissipation  in  case  S2  to  —  8.2,  bringing  it  much  closer  to  the  value  obtained  in  case  E2. 

The  resolved  and  SGS  shear  stresses  are  plotted  in  figure  15.  The  SGS  shear  stresses  for  the  Smagorinsky  model 
and  the  dynamic  model  are  concentrated  in  the  wall  regions  and  are  close  to  zero  outside  of  the  wall  boundary  layers. 
The  estimation  procedure  gives  SGS  stress  distributed  with  a  fairly  constant  slope  across  the  entire  channel  width 
and  a  required  decay  to  zero  at  the  walls.  The  behavior  of  all  models  in  this  case  is  consistent  with  the  behavior 
observed  in  the  low  Reynolds  number  cases  (see  Figs.  1(d)  and  6(d)),  which  favored  predictions  of  the  estimation 
procedure.  The  large  values  of  the  SGS  stress  concentrated  in  the  wall  region  for  the  Smagorinsky  model  and  the 
dynamic  model  are  the  result  of  the  dependence  of  the  eddy  viscosity  formulation  on  the  resolved  rate-of-stress  tensor 
which  is  dominated  by  large  vertical  gradients  of  the  streamwise  velocity  in  the  boundary  layers.  Adjusting  the 
Smagorinsky  constant  from  0.1  to  0.06  would  lower  the  peak  values  in  case  S2  from  0.6  to  0.22  which  is  much  closer  to 
the  dynamic  model  value.  The  differences  in  model  predictions  for  the  SGS  shear  stress  lead  to  substantial  differences 
in  predicted  profiles  of  the  resolved  shear  stress  shown  in  Fig.  15(b).  This  is  because  in  a  steady  state  the  total  shear 
stress  (sum  of  viscous,  resolved,  and  SGS  components)  depends  linearly  on  the  wall  normal  coordinate  z.  Since  the 
mean  velocities  are  similar  in  all  cases,  the  viscous  stresses  are  also  similar  and  differences  in  the  modeled  SGS  stresses 
must  be  compensated  by  differences  in  the  resolved  stresses. 

The  results  from  the  high  Reynolds  number  LES  cases  confirm  conclusions  about  the  estimation  procedure  drawn 
previously  from  a  priori  analysis  and  LES  for  low  Reynolds  number  turbulence.  The  procedure  predicts  mean  and 
fluctuating  velocities  in  general  agreement  with  standard  eddy  viscosity  models,  gives  better  representation  of  SGS 
stresses,  and  ew:counts  naturally  for  backscatter  without  any  numerical  instabilities. 


V.  CONCLUSIONS 

In  this  paper,  we  have  proposed  a  new  approach  to  the  problem  of  subgrid-scale  modeling  where  one  attempts  to 
find  an  estimation  of  a  velocity  field  which  contains  a  limited  range  of  subgrid-scales  unknown  in  LES  and  which 
is  consistent  with  the  resolved  velocity.  The  estimated  velocity  is  then  used  to  compute  the  SGS  stress  from  the 
definition.  Thus,  in  the  proposed  method  focus  is  not  on  devising  an  expression  for  the  SGS  stress  tensor,  but  on 
modeling  unknown  subgrid-scales.  Using  results  from  previous  a  priori  analyses  of  DNS  databases,  it  was  concluded 
that  only  subgrid-scales  which  are  at  most  a  factor  of  two  smaller  than  the  smallest  resolved  scales  need  to  be 
modeled.  This  distinguishes  the  proposed  method  from  the  traditional  eddy  viscosity  concepts  which  rely  implicitly 
on  the  assumption  of  wide  separation  between  interacting  scales.  The  subgrid-scales  are  determined  using  kinematic 
compatibility  conditions  between  the  unknown  and  the  resolved  field  and  a  dynamic  condition  reflecting  the  influence 
of  nonlinear  interactions  between  resolved  scales  on  subgrid-scales.  Through  the  kinematic  condition  the  estimated 
velocity  and  in  turn  the  SGS  stress  tensor  depend  explicitly  on  the  filtering  procedure  used.  This  is  an  attractive 
feature  of  the  method,  because  the  SGS  stress  computed  from  a  given  fully  resolved  velocity  field  must  depend  on  the 
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filter  and  many  SGS  models  fail  to  exhibit  such  a  dependence.  For  example,  the  classical  Smagorinsky  expression  is 
used  in  exactly  the  ssme  form  for  the  LES  equations  obtained  with  different  filters.  The  velocity  estimation  procedure 
has  been  implemented  in  the  current  study  only  for  the  tophat  filter  and  for  fields  which  can  be  represented  in 
terms  of  Fourier  expansions.  Since  the  explicit  properties  of  the  tophat  filter  and  of  the  discrete  Fourier  transforms 
were  used,  the  method  in  its  current  form  cannot  be  applied  directly  for  other  combinations  of  the  filter  and  the 
numerical  representation  of  the  velocity  field.  The  expression  for  the  SGS  stress  does  not  depend  on  any  adjustable 
constants,  because  beyond  the  need  for  an  additional  filter  width,  no  adjustable  constants  are  introduced  in  the 
velocity  estimation  procedure.  By  construction,  the  SGS  stress  is  Galilean  invariant  and  has  the  correct  near  wall 
behavior. 

The  proposed  SGS  model  was  implemented  in  a  channel  flow  code  and  evaluated  in  both  a  priori  analysis  and 
actual  large  eddy  simulations  at  two  different  Reynolds  numbers.  The  results  were  compared  with  results  provided 
by  LES  performed  with  the  Smagorinsky  model  and,  at  higher  Reynolds  number,  with  the  dynamic  model  results  of 
Piomelli  [51].  The  estimation  procedure  was  found  to  predict  mean  and  rms  velocities  with  accuracy  comparable  to 
the  eddy  viscosity  models.  However,  the  predictions  for  the  SGS  stresses  themselves  were  in  much  better  qualitative 
and  quantitative  agreement  with  the  exsct  results  than  those  obtained  with  the  Smagorinsky  model.  The  non-eddy 
viscosity  models  usually  encounter  problems  with  inadequate  SGS  dissipation,  forcing  the  addition  of  dissipative  eddy 
viscosity  expressions  in  actual  LES.  The  estimation  procedure  in  low  Reynolds  number  LES  was  found  to  produce 
SGS  dissipation  comparable  to  the  Smagorinsky  model.  From  the  results  for  the  high  Reynolds  number  case,  we  can 
infer  that  the  estimation  procedure  gives  SGS  dissipation  comparable  to  the  dissipation  for  the  dynamic  model  and 
the  Smagorinsky  model  with  a  constant  Cs  «  0.06.  In  all  cases  the  net  SGS  dissipation  is  sufficient  to  guarantee  the 
correct  evolution  of  the  flow.  The  new  model  also  predicts  significant  local  backscatter  without  any  adverse  effects  on 
the  numerical  stability  of  the  simulations.  The  above  results  strongly  suggest  that  the  estimation  procedure  proposed 
and  investigated  in  this  paper  is  a  physically  and  numerically  viable  approach  to  the  problem  of  SGS  modeling. 

The  proposed  modeling  procedure  should  be  investigated  further  to  answer  several  important  questions.  The 
procedure  is  not  unique,  because  there  are  many  possible  estimated  fields  which  will  be  consistent  with  a  resolved 
field.  What  is  the  optimal  procedure  for  choosing  from  such  possible  fields?  The  simulations  with  the  estimation  model 
at  higher  Reynolds  number  were  started  with  the  initial  condition  taken  from  the  lower  Reynolds  number  case,  and 
then  run  until  steady  state  was  reached.  However,  the  transient  evolution  was  not  studied.  How  will  the  model  behave 
in  non-stationary  and  nonequilbrium  flows?  The  estimated  velocity  is  represented  on  a  mesh  with  twice  as  many  grid 
points  as  the  resolved  velocity  in  each  Cartesian  direction  in  which  filtering  is  applied.  Such  a  representation  could 
be  numerically  expensive  in  three  dimensional  simulations  if  filtering  is  applied  in  all  three  directions.  Is  it  possible 
to  obtain  good  results  with  fewer  points?  Can  the  procedure  be  sequentially  applied  in  each  direction  rather  than 
simultaneously  in  all  directions  as  done  in  this  work?  For  very  high  Reynolds  numbers  nonlocal  interactions  between 
the  resolved  scales  and  subgrid-scales  in  the  range  k  >  2kc  no  longer  will  be  negligible.  Can  their  effect  be  properly 
represented  by  a  simple  eddy  viscosity  because  of  the  separation  of  scales?  Currently,  LES  with  the  estimation  model 
require  about  40%  more  time  to  run  than  LES  with  the  standard  Smagorinsky  model.  Can  this  timing  be  improved? 
Finally,  an  important  question  is  how  to  implement  the  modeling  principles  in  general  finite  difference  codes.  We  plan 
to  address  the  above  questions  in  future  work. 
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FIG.  1.  The  plane  averaged  components  of  the  SGS  stress  tensor  modified  by  subtracting  ^SijTkk-  Solid  line:  DNS;  dashed 
line:  the  Smagorinsky  model;  dashed-dotted  line:  the  velocity  estimation  model,  (a)  <  tu  —  -t**  >;  (b)  <  n2  —  ^Tkk  >;  (c) 
<  ^33  -  ^Tkk  >;  (d)  <  ri3  >. 


FIG.  2.  The  plane  averaged  components  of  the  SGS  force  (41).  Solid  line:  DNS;  dashed  line:  the  Smagorinsky  model; 
dashed-dotted  line:  the  velocity  estimation  model,  (a)  <  ^(rij  -  jSijTkk)  >;  (b)  <  -^{nj  —  ^SsjTkk)  >. 


FIG.  3.  SGS  dissipation.  Solid  line:  DNS;  dashed  line:  the  Smagorinsky  model;  dashed-dotted  line:  the  velocity  estimation 
model,  (a)  The  local  SGS  dissipation  averaged  over  horizontal  planes,  (b)  Negative  values  of  the  local  SGS  dissipation  (forward 
transfer)  and  positive  values  (bcickscatter)  averaged  over  horizontal  planes. 


FIG.  4.  The  plane  correlation  coefficients  between  the  DNS  and  modeled  components  of  T13.  Solid  line:  the  Smagorinsky 
model;  dashed-dotted  line:  the  velocity  estimation  model. 


FIG.  5.  The  plane  correlation  coefficients  between  the  SGS  stress  component  T13  and  the  resolved  rate-of-strain  tensor 
component  5 13.  Solid  line:  DNS;  dashed  line:  quantities  modeled  using  the  Smagorinsky  model;  dashed-dotted  line:  quantities 
modeled  using  the  velocity  estimation  model. 


FIG.  6.  The  SGS  stress  tensor  modified  by  subtracting  ^SijTkk-  Solid  line:  DNS;  dashed  line:  case  Si;  dashed-dotted-dotted 
line:  case  El.  (a)  <  m  -  ^Tkk  >;  (b)  <  T22  -  ^Tkk  >;  (c)  <  T33  -  jTkk  >;  (d)  <  m  >. 


FIG.  7.  The  plane  averaged  components  of  the  SGS  force.  Solid  line:  DNS;  dashed  line:  case  SI;  dashed-dotted-dotted  line: 
case  El.  (a)  <  ^(nj  -  ^Sijnk)  >;  (b)  <  -  ^SsjTkk)  >■ 
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FIG.  8.  SGS  dissipation.  Solid  line:  DNS;  dashed  line:  case  SI;  dashed-dotted-dotted  line:  case  El.  (a)  The  loceii  SGS 
dissipation  averaged  over  horizontal  planes,  (b)  Negative  values  of  the  local  SGS  dissipation  (forward  transfer)  and  positive 
values  (backscatter)  averaged  over  horizontal  planes. 

FIG.  9.  The  plane  correlation  coefficients  between  the  SGS  stress  component  ns  and  the  resolved  rate-of-strain  tensor 
component  5 13.  Solid  line:  DNS;  dashed  line:  case  Si;  dashed-dotted-dotted  line:  case  El. 

FIG.  10.  Mean  velocity  profiles.  Solid  line:  DNS;  dashed  line:  case  SI;  dashed-dotted-dotted  line:  case  El;  dotted  line:  case 
A;  dashed-dotted  line:  =  2.51n(z‘*’)  +  5.5;  long  dashed  line:  =  z"*". 

FIG.  11.  Rms  turbulent  velocity.  Solid  line:  DNS  (filtered);  dashed-dotted  line:  DNS  (unfiltered);  dashed  line:  case  Si; 
dashed-dotted-dotted  line:  case  El;  dotted  line:  case  A.  (a)  Streamwise  component  Urm«.  (b)  Spanwise  component  Vrma-  (c) 
WaU-normal  component  Wrma- 

FIG.  12.  Mean  velocity  profiles.  Solid  line:  Piomelli  (1993);  dashed  line:  S2;  dashed-dotted-dotted  line:  E2;  dashed-dotted 
line:  u"*"  =  2.51n(z'*')  -|-  5.5;  long  dashed  line:  =  z'*’. 

FIG.  13.  Rms  turbulent  velocity.  Solid  line:  Piomelli  (1993);  dashed  line:  S2;  dashed-dotted-dotted  line:  E2.  (a)  Streamwise 
component  Wrms.  (b)  Sp£uiwise  component  Vrms-  (c)  Wall-normcil  component  Wrms. 

FIG.  14.  The  decomposition  of  SGS  dissipation  averaged  over  horizontal  planes  into  forward  transfer  (negative  values)  and 
backscatter  (positive  values).  Dashed  line:  case  S2;  dashed-dotted-dotted  line:  case  E2.  The  pezJc  in  the  forward  transfer  for 
case  S2  is  fa  —256. 

FIG.  15.  Shear  stress  profiles.  Solid  line:  Piomelli  (1993);  dashed  line:  case  S2;  dashed-dotted-dotted  line:  case  E2.  (a)  SGS 
stress  component  <  ri3  >.  (b)  Resolved  stress  component  <uw>. 
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TABLE  1.  The  global  correlation  coefficients  between  exact  and  modeled  SGS  quantities.  EM-estimation  model; 
EMl-estimation  model  with  modified  stresses  Tij  —  SM-Smagorinsky  model. 


Quantity 

EM 

EMI 

SM 

Til 

0.53 

0.52 

0.12 

T12 

0.44 

0.44 

0.07 

T13 

0.51 

0.51 

0.31 

T22 

0.54 

0.46 

-0.003 

T23 

0.32 

0.32 

0.06 

T33 

0.52 

0.55 

0.16 

NfGS 

0.15 

0.14 

0.13 

0.12 

-0.02 

-0.02 

0.25 

0.43 

0.12 

rpSGS 

0.15 

0.14 

0.1 

esGS 

0.4 

0.4 

0.38 

TABLE  II.  The  parameters  implemented  in  the  numerical  simulations 


Case 

Grid 

Rer 

time 

Model 

A 

32  X  32  X  65 

180 

4.2 

No  model 

SI 

32  X  32  X  65 

180 

4.2 

Smagorinsky 

El 

32  X  32  X  65 

180 

4.2 

Estimation 

S2 

48  X  64  X  65 

1050 

6.2 

Smagorinsky 

E2 

48  X  64  X  65 

1050 

22.2 

Estimation 

TABLE  III.  The  integrated  values  of  SGS  dissipation. 


Case 

Rct 

<^SGS> 

net 

210 

-8.3 

0.7 

-7.6 

SI 

180 

-7.1 

0.0 

-7.1 

El 

180 

-6.8 

1.8 

-5.0 

S2 

1050 

-22.7 

0.0 

-22.7 

E2 

1050 

-16.9 

6.3 

-10.6 
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FIG.  1.  The  plame  averaged  components  of  the  SGS  stress  tensor  modified  by  subtracting  ^SijTkk-  Solid  line:  DNS;  dashed 
line:  the  Smagorinsky  model;  dashed-dotted  line:  the  velocity  estimation  model,  (a)  <  rn  —  ^Zkk  >;  (b)  <  T22  —  ^Tkk  >;  (c) 
<  T33  ~  ^Tkk  >;  (d)  <  ri3  >. 


DNS 

Smagorinsky  model 
Velocity  estimation  model 


FIG.  2.  The  plane  averaged  components  of  the  SGS  force  (41).  Solid  line:  DNS;  dashed  line:  the  Smagorinsky  model 
dashed-dotted  line:  the  velocity  estimation  model,  (a)  <  -^{nj  —  >;  (b)  <  >. 


FIG.  3.  SGS  dissipation.  Solid  line:  DNS;  dashed  line:  the  Smagorinsky  model;  dashed-dotted  line:  the  velocity  estimation 
model,  (a)  The  local  SGS  dissipation  averaged  over  horizontal  planes,  (b)  Negative  VcJues  of  the  local  SGS  dissipation  (forward 
transfer)  and  positive  values  (bau:kscatter)  averaged  over  horizontal  planes. 
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Velocity  estimation  model 
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FIG.  5.  The  plane  correlation  coefficients  between  the  SGS  stress  component  ris  and  the  resolved  rate-of-strain  tensor 
component  S13.  Solid  line:  DNS;  dashed  line:  quantities  modeled  using  the  Smagorinsky  model;  dashed-dotted  line:  quantities 
modeled  using  the  velocity  estimation  model. 


FIG.  6.  The  SGS  stress  tensor  modified  by  subtracting  jSijTkk-  Solid  line:  DNS;  dctshed  line:  case  SI;  dashed-dotted-dotted 
line:  case  El.  (a)  <  tu  -  ^Tkk  >;  (b)  <  m  -  ^Tkk  >;  (c)  <  T33  -  ^Tkk  >;  (d)  <  ns  >. 
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FIG.  6(c) 
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FIG.  6(d) 
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FIG.  7.  The  plane  averaged  components  of  the  SGS  force.  Solid  line:  DNS;  dashed  line:  case  SI;  dashed-dotted-dotted  line: 
case  El.  (a)  <  -  iJijTkfc)  >;  (b)  <  -  |<53,T*t)  >. 


33 


FIG.  11.  Rms  turbulent  velocity.  Solid  line:  DNS  (filtered);  dashed-dotted  line:  DNS  (unfiltered);  dashed  line:  case  SI; 
dashed-dotted-dotted  line:  case  El;  dotted  line:  case  A.  (a)  Streeimwise  component  Urma-  (b)  Spanwise  component  Vrms.  (c) 
Wall-normcJ  component  Wrma- 


FIG.  15(b) 
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Abstract 

Using  previously  established  properties  of  nonlinear  subgrid-scale  interactions  we  derive 
deterministic  models  of  backscatter  for  large  eddy  simulations.  The  models  affect  only  weakly 
the  total  kinetic  energy  of  turbulence  providing  approximately  equal  amounts  of  forward  and 
back-scatter.  The  modeling  procedures  are  tested  on  a  simple  case  of  forced  isotropic  turbu¬ 
lence  with  the  inertial  range  and  in  turbulent  channel  flow. 
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1  Introduction 


Currently  several  subgrid-scale  (SGS)  models  are  commonly  used  in  large  eddy  simulations  (LES) 
of  turbulence.  The  Smagorinsky  model  ^  and  the  models  for  homogeneous  turbulence  proposed  by 
Kraichnan  ^  and  Chollet  and  Lesieur  ^  have  a  status  of  classical  models.  These  models  are  also 
used  as  a  starting  point  to  construct  newer  models:  the  dynamic  model  of  Germano  et  al.  ^  and  its 
variations  and  the  structure  function  model  of  Metals  and  Lesieur  A  recent  review  of  these 
as  well  as  several  other  SGS  models  is  given  by  Lesieur  The  common  and  important  feature  of 
these  successful  SGS  models  is  that  they  properly  account  for  the  global  subgrid-scale  dissipation, 
i.e.  the  net  energy  flux  from  the  resolved  to  the  subgrid  scales. 


The  models  are  less  successful  in  accounting  for,  or  often  do  not  attempt  to  account  directly 
for  the  phenomenon  of  the  inverse  energy  transfer  from  the  subgrid  scales  to  the  resolved  scales. 
This  process,  often  called  backscatter,  is  implicitly  interpreted  in  two  different  ways.  In  the  context 
of  the  analytical  theories  of  turbulence  Leslie  and  Quarini  noted  that  the  resolved  scale  energy 
equation  contains  two  distinct  terms,  one  describing  energy  drain  from  the  resolved  scales  to  the 
subgrid  scales  and  the  other,  the  backscatter  term,  describing  energy  transfer  from  the  subgrid 
scales  to  the  resolved  scales.  The  total  energy  transfer  is  the  sum  of  these  two  terms  and  when  it  is 
represented  in  a  form  of  the  eddy  viscosity  the  resulting  net  eddy  viscosity  contains  contributions 
from  both  forward  and  inverse  energy  transfers.  Therefore  such  a  net  eddy  viscosity  takes  indirectly 
into  account  the  presence  of  the  backscatter  by  decreasing  the  value  of  the  eddy  viscosity  associated 
with  the  forward  energy  transfer  term.  A  somewhat  different  meaning  is  given  to  the  backscatter 
term  in  the  context  of  a  priori  analyses  of  direct  numerical  simulations  data.  In  isotropic  turbulence 
simulations  the  rate  of  energy  change  of  an  individual  wavenumber  mode  k  may  be  calculated  and 
is  either  negative  or  positive.  The  latter  case  corresponds  to  backscatter  local  in  a  wavenumber. 
When  negative  and  positive  contributions  are  averaged  over  wavenumber  shells  of  radius  k  the 
total  transfer  is  split  into  forward-  and  back-scatter,  respectively,  formally  corresponding  to  the 
decomposition  used  in  the  analytical  theories  of  turbulence.  However,  it  has  not  been  shown  yet 
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that  the  backscatter  calculated  in  DNS  is  indeed  approximated  by  the  backscatter  term  in  the 
analyticcd  theories.  In  fact  the  spectrum  of  backscatter  computed  from  low  Reynolds  number 
turbulence  simulations  by  Domaxadzki  et  al.  does  not  have  properties  expected  from  theoretical 
analysis.  Therefore  caution  must  be  used  in  identifying  backscatter  in  DNS  with  the  backscatter 
term  in  the  analytical  theories  until  detailed  comparison  for  different  flows  is  made.  Another 
important  distinction  is  that  the  theoretical  net  eddy  viscosity  is  generally  positive  implying  that  it 
decreases  the  amplitudes  of  individual  modes  k  while  simulations  show  that  individual  modes  may 
grow  or  decay  because  of  the  subgrid-scale  interactions.  Similarly,  in  a  physical  space  representation 
the  actual  subgrid-scale  dissipation  is  locally  frequently  negative  while  the  classical  eddy  viscosity 
models  always  give  positive  values.  It  also  must  be  noted  that  there  is  no  simple  relation  between 
local  backscatter  in  physical  and  spectral  representation.  Therefore,  the  backscatter  term  is  used 
to  qualitatively  describe  in  all  cases  the  same  phenomenon,  energy  transfer  from  the  subgrid  to 
the  resolved  scales,  but  one  has  to  be  aware  that  in  all  of  those  case  the  quantitative  definitions  of 
backscatter  are  not  identical. 

Analyses  of  direct  numerical  simulations  of  various  turbulent  flows  1142,13,14  established  that 
the  backscatter  is  comparable  to  and  often  larger  than  the  net  subgrid-scale  dissipation.  This 
clearly  implies  that  accounting  for  backscatter  is  desirable  in  LES  on  purely  physical  grounds. 
However  the  motivation  for  modeling  backscatter  phenomena  in  typical  applied  turbulence  LES 
is  lacking.  This  is  because  important  mean  and  rms  turbulent  quantities  appear  to  be  sensitive 
only  to  the  net  SGS  dissipation.  However,  as  argued  by  Schumann  it  is  likely  that  backscatter 
effects  become  more  important  in  situations  where  subgrid  scales  contain  a  significant  fraction  of 
the  total  energy,  or  more  generally,  when  the  resolution  is  necessarily  coarse  as  in  all  geophysical 
applications.  Despite  the  lack  of  a  strong  practical  motivation  at  the  present  time  for  modeling 
backscatter,  its  presence  in  turbulent  flows  constitutes  a  sufficient  physical  reason  that  motivates 
work  in  this  area.  Backscatter  can  be  introduced  into  LES  by  using  random  forcing  (Leith 
Chasnov  Schumann  or  deterministically  by  extrapolating  from  the  dynamics  of  the  resolved 
scales  in  the  framework  of  the  dynamic  models  as  discussed  by  Carati  et  al.  or  in  the  framework 
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of  the  similaxity  models  stochastic  models  can  be  designed  to  provide  proper  amounts 

of  net  inverse  energy  transfer  and  the  fact  that  the  subgrid  scales  are  not  known  makes  random 
forcing  a  natural  choice.  However,  as  observed  by  Schumann  the  subgrid  scales  being  solutions 
of  Navier-Stokes  equations  can  not  be  completely  random.  The  purely  stochastic  models  are  also 
unable  to  describe  in  a  physically  satisfactory  manner  the  nonlinear  interactions  between  highly 
correlated  (smallest  resolved  and  largest  subgrid)  scales  which  are  responsible  for  the  subgrid-scale 
energy  transfer  The  effects  of  such  strong  correlations  axe  observed  by  Haxtel  and  Kleiser 

in  the  wall  region  of  channel  flow  where  the  subgrid-scales  may  constitute  with  the  resolved  scales 
part  of  the  same  physical  coherent  structure.  Moreover,  to  be  consistent,  if  the  effect  of  the  subgrid 
scales  on  the  resolved  scales  is  entirely  random  then  not  only  backscatter  but  also  forwaxd  transfer 
should  be  modeled  in  a  non-deterministic  fashion,  and  LES  practice  indicates  that  it  is  not  necessary. 
Therefore,  while  some  level  of  stochasticity  in  a  SGS  model  is  desirable,  deterministic  models  offer 
a  reasonable  alternative  in  view  of  a  lack  of  detailed  information  about  statistical  properties  of  SGS 
interactions. 

1  O 

The  backscatter  in  the  dynamic  model  has  a  good  physical  interpretation  but  its  origin  in  the 
eddy  viscosity  model  appears  to  limit  backscatter  to  values  lower  than  observed  in  a  ‘priori  tests. 
The  similarity  models  provide  the  most  natural  way  of  introducing  backscatter  deterministically. 
However,  in  this  approach  one  usually  attempts  to  model  the  entire  subgrid-scale  stress,  possibly  a 
more  difficult  task  than  modeling  backscatter  separately  and  then  combining  it  with  a  dissipative 
expression  (giving  a  version  of  a  mixed  model).  Moreover,  such  models  often  use  formal  similarities 
between  subgrid  stresses  computed  with  different  filter  widths  rather  than  known  properties  of  the 
subgrid-scale  interactions.  As  a  result,  for  instance,  the  SGS  stress  tensor  in  the  Bardina  model 
vanishes  if  the  sharp  spectral  filter  is  employed.  For  other  filters  the  stress  tensor  is  nonzero  but  it 
significantly  underestimates  the  real  energy  transfers  as  reported  by  Liu  et  al. 

These  difficulties  and  uncertainties  in  accounting  for  backscatter  in  LES  suggest  that  better 
models  may  be  expected  only  if  more  information  about  physics  of  this  process  is  employed  in  the 
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modeling.  The  goal  of  this  paper  is  to  propose  procedures  for  modeling  backscatter  effects  which  are 
based  on  properties  of  nonlinear  subgrid  scale  interactions  observed  in  direct  numerical  simulations 
of  turbulence. 

2  Backscatter  models 

Using  a  priori  analysis  of  isotropic  turbulence  data,  Kerr  et  al.  have  shown  that  most  of  the 
backscatter  effects  are  described  by  the  term  obtained  from  the  decomposition  of  the  full  nonlinear 
term  in  the  vorticity  form: 


B  =  (u^  X  (1) 

where  u  and  u  are  the  velocity  and  the  vorticity  fields  and  the  superscript  C  signifies  filtering  to 
large  scales,  resolved  in  LES,  and  S  signifies  filtering  to  subgrid  scales,  unknown  in  LES.  If  a  tilde 
is  used  to  denote  filtering  to  the  large  scales  the  above  formula  becomes 


B  =  u  X  (w  —  (D).  (2) 

More  specifically  Kerr  et  al.  found  that  term  (1)  provides  large  but  about  equal  amounts 
of  forward  and  back-scatter  such  that  it  contributes  very  little  to  the  net  SGS  energy  transfer. 
The  remaining  terms  in  the  decomposition  are  responsible  for  the  majority  of  the  net  transfer. 
These  results  suggest  that  the  effects  of  backscatter  can  be  included  in  LES  by  adding  to  a  purely 
dissipative  SGS  model  which  properly  accounts  for  the  net  transfer,  an  additional  term  which  is 
based  on  the  above  expression.  This  is  the  strategy  that  we  explore  in  this  paper. 

Formulas  (1)  and  (2)  give  the  vector  product  of  the  resolved  velocity  with  the  unresolved  vor¬ 
ticity,  the  result  filtered  subsequently  to  the  resolved  scales.  To  model  this  expression  for  the  use  in 
LES  we  need  a  model  for  the  subgrid  vorticity.  Two  such  models  based  on  the  small-scale  vorticity 
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production  axe  discussed  by  Kerr  et  al.  but  we  consider  another  approach  here.  The  simplest 
model  is  based  on  the  similarity  concept.  The  resolved  scales  are  divided  into  large  and  small  scales 
by  an  appropriate  filtering  operation  and  the  above  formulas  are  applied.  If  this  filtering  operation, 
often  called  “a  test  filtering”  in  the  dynamic  SGS  models,  is  denoted  by  the  overbar,  the  model  is 


B(^)  =  ux(u;-u;).  (3) 

In  what  follows  we  will  simplify  the  notation  by  dropping  tilde,  i.e.  assuming  implicitly  that 
quantities  without  a  tilde  axe  the  resolved  LES  quantities  and  the  superscript  C  and  the  overbax 
signify  the  test  filtering.  The  model  (3)  is  then  written  simply  as  =  u  x  (a>  —  W). 

Generally,  a  subgrid-scale  model  in  a  form  of  the  divergence  of  a  second-order  tensor  is  preferred 
because  this  form  assures  that  the  mean  momentum  of  the  flow  is  affected  by  the  subgrid-scale 
interactions  only  at  the  boimdaries  as  required  by  the  form  of  the  SGS  stress  tensor  in  the  LES 
equations.  In  particular,  subgrid-scale  interactions  cannot  generate  a  mean  flow  in  isotropic  turbu¬ 
lence.  This  condition  is  satisfied  by  (3)  only  if  the  sharp  spectral  filter  is  used  as  was  the  case  in 
the  a  priori  analysis  of  Kerr  et  al.  If  other  filtering  operations  are  to  be  used,  the  model  (3) 
should  be  replaced  by  a  model  in  the  stress  tensor  form.  The  expression  (1)  can  be  written  as 


Bi  = 


^  /  S  C\C  ,  /  C  ^  S\C 


(4) 


where  the  first  term  has  a  desired  form  suggesting  using  it  as  a  model.  This  implies  the  backscatter 


stress  tensor  in  a  form  Tjj  =  (ufuf)^  =  (u,-  —  Ui)uj 


The  above  models  originated  in  the  analyses  of  homogeneous,  isotropic  turbulence  where  the 
mean  velocity  vanishes  and  u(x,<)  is  the  fluctuating  velocity.  If  the  mean  flow  is  non-zero  and 
is  included  in  the  definition  of  u(x,  t)  then  both  models  will  not  be  Galilean  invariant  because  of 
proportionality  to  the  large  scale  velocity  uf.  One  way  to  recover  the  Galilean  invariance  of  the 
models  for  cases  with  nonvanishing  mean  flow  and  to  maintain  their  consistency  with  the  original 
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analyses  is  to  apply  the  modeling  expressions  to  the  fluctuating  components  of  the  velocity  only. 
This  approach  is  favored  for  the  model  (3)  applied  to  homogeneous  flows.  For  inhomogeneous  flows 
the  divergence  form  is  preferred  and  in  that  case  another  standard  method  to  maJce  the  model 
Galilean  invariant  is  to  redefine  SGS  stresses  as  proposed  by  Germano 


T^p  =  ufuf-  uf  uf  =  (Uf  -  Ui)Uj  -  (Ui  -  Ui)  Uj 
with  the  corresponding  model 


= 


..(2) 


(5) 


(6) 


Physically  models  (3)  and  (6)  axjcount  for  modifications  of  the  large  resolved  scales  C  by  the  small 
resolved  scales  S  while  the  small  scales  are  being  advected  by  the  large  ones. 


In  principle  both  formulas  could  also  contain  a  multiplicative  constant  and  another  free  param¬ 
eter  is  the  ratio  of  the  filter  width  used  in  formulas  (3)  and  (6)  to  the  filter  width  of  the  LES.  In 
the  original  backscatter  expression  (1),  the  multiplicative  constant  is  equal  to  unity  and  we  choose 
to  keep  this  value  as  well  in  the  modeling  expressions  (3)  and  (6).  Therefore,  the  filter  width 
ratio  remains  as  the  only  free  parameter.  The  models  affect  only  the  largest  resolved  scales  deter¬ 
mined  by  the  filter  C  while  it  is  known  that  the  backscatter  interactions  affect  all  resolved  scales  in 
LES  'p^g  suggests  that  the  range  of  scales  S  should  be  made  as  small  as  possible  while 

using  these  models.  An  analysis  of  DNS  data  for  low  Reynolds  number  flows  and  the  results  of  LES 
reported  in  the  next  section  indicate  that  the  test  filter  cutoff  wavenumber  should  be  in  the  range 
between  50  and  90%  of  the  LES  cutoflf  value. 


3  Large  eddy  simulations  with  backscatter  models 


Analyses  of  the  direct  numerical  simulations  data  by  themselves  only  suggest  but  do  not  guarantee 
that  the  desirable  features  of  the  models  will  be  in  fact  observed  in  actual  large  eddy  simulations. 
In  this  section  we  assess  the  backscatter  models  explicitly. 
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3.1  Forced  isotropic  turbulence 


A  simple  case  of  isotropic  turbulence  with  the  inertial  subrange  was  chosen  to  see  if  the  models 
provide  approximately  equal  amounts  of  forward  and  inverse  transfers  in  LES  without  affecting 
significantly  the  net  SGS  transfer.  First,  we  performed  pseudo-spectral  large  eddy  simulations  of 
forced  isotropic  turbulence  on  a  32®  grid.  We  employed  the  forcing  scheme  of  Sullivan  et  al.  which 
keeps  the  total  energy  of  several  low  wavenumber  modes  constant,  but  allows  evolution  of  individual 
modes  through  nonlinear  interactions,  subject  to  the  global  energy  constraint.  The  subgrid-scale 
model  used  was  a  constant  spectral  eddy  viscosity  discussed  by  Lesieur  ^ 


i/T  =  0A02y/Eikm)lkm,  (7) 

where  E{k)  is  the  energy  spectrum  and  km  is  the  maximum  wavenumber  in  LES,  i.e.  the  wavenum¬ 
ber  marking  the  boundary  between  resolved  and  subgrid  scales.  The  constant  0.402  in  (7)  is  obtained 
assuming  the  inertial  range  with  the  Kolmogoroff  constant  Ko  =  1.4.  The  program  was  run  until 
statistically  stationary  state  was  reached.  The  molecular  viscosity  term  was  retained  in  the  simula¬ 
tions  and  in  the  stationary  state  j/t/i/  w  500.  After  a  stationary  state  was  reached,  the  simulations 
were  run  for  additional  several  hundred  time  steps  keeping  the  eddy  viscosity  model  (7)  and  adding 
backscatter  models  (3)  and  (6)  to  the  LES  equations. 

In  Fig.  1  we  plot  the  spectral  Kolmogoroff  function  defined  as 

Ko{k)  =  E{k)l{^l^k-^l^)  (8) 

for  the  base  LES  with  the  eddy  viscosity  model  and  for  runs  with  the  additional  backscatter  models. 
In  (8),  e  is  the  energy  flux  down  the  spectrum,  equal  to  the  energy  input  into  the  system  necessary 
to  maintain  the  steady  state.  For  the  perfect  inertial  range  spectrum  consistent  with  the  eddy 
viscosity  model  used,  (8)  should  be  constant  equal  to  1.4.  Ko{k)  in  the  simulations  varies  between 
1  and  2.2  showing  similar  trends  for  all  models  used  though  adding  the  backscatter  models  results 
in  slightly  improved  predictions.  The  level  of  variation  in  Ko  is  similar  to  that  observed  in  typical 
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LES  of  isotropic  turbulence  and  much  better  results  should  not  be  expected  in  view  of  the  fairly 
simple  eddy  viscosity  model  used.  The  important  feature  observed  in  Fig.  1  is  that  all  backscatter 
models  used  do  not  significantly  affect  the  energy  spectra  providing  results  comparable  to  those 
obtained  using  the  eddy  viscosity  model.  However,  the  eddy  viscosity  has  a  purely  dissipative  effect 
on  the  energy  of  the  flow  resulting  in  the  energy  dissipation  rate  equal  to  the  the  energy  input  rate  e 
in  the  steady  state.  The  role  of  the  backscatter  terms  is  to  model  a  (positive)  influx  of  energy  from 
the  subgrid  to  the  resolved  scales.  In  Fig.  2  we  plot  as  a  function  of  time  (positive)  backscatter  and 
(negative)  forward  transfer  components  for  the  energy  flux  to/from  the  resolved  scales  associated 
with  each  model  and  normalized  by  the  global  transfer  rate  e.  These  energy  transfer  rates  were 
computed  by  summing  respective  positive  and  negative  values  of  contributions  that  the  models 
make  to  the  kinetic  energy  equation  in  the  physical  space.  In  all  cases  the  test  filtering  was  done 
with  the  sharp  spectral  filter  in  each  direction,  x,  y,  and  z,  with  the  cutoff  wavenumber  kc  =  ZlAkm- 
Both  models  provide  approximately  equal  amounts  of  forward-  and  back-scatter,  between  60%  and 
80%  of  the  global  energy  dissipation  rate,  depending  on  the  model.  These  values  compare  favorably 
with  values  observed  in  a  priori  analyses  of  DNS  data,  but  are  significantly  larger  than  found  by 
Carati  et  al.  for  the  current  state-of-the-art  dynamic  DLM(k)  model. 

3.2  Channel  flow  turbulence 

We  have  performed  large  eddy  simulations  of  turbulent  channel  flow  using  a  pseudo-spectral  Fourier- 
Chebyshev  numerical  method  with  the  resolution  of  32  x  32  modes  in  the  streamwise  x  and  the 
spanwise  y  directions  and  65  modes  in  the  vertical  z  direction.  The  physical  parameters  in  the  sim¬ 
ulations  were  chosen  to  match  parameters  used  in  direct  numerical  simulations  of  the  same  problem 
performed  by  Kim  et  al.  at  Reynolds  number  Re  =  3300  based  on  the  mean  centerline  velocity 
and  the  channel  half  width  h.  The  corresponding  Reynolds  number  based  on  friction  velocity  Ur  and 
h  is  approximately  Rcr  =  180.  The  LES  results  were  compared  with  a  DNS  database  of  Gilbert 
made  availabe  to  us  by  Professor  L.  Kleiser  (see  also  Gilbert  and  Kleiser  and  Hartel  et  al.  ^  ).  In 
these  DNS  channel  flow  was  simulated  with  the  resolution  of  160^(horizontal)xl29(vertical)  modes 
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Kolmogoroff  Constant 


Figure  1:  Predictions  of  the  Kolmogoroff  function.  Solid  line:  eddy  viscosity  (7)  without  backscatter; 
broken  line:  backscatter  model  (3);  dotted  line:  backscatter  model  (6). 
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at  slightly  higher  Reynolds  number  Rcr  =  210.  Turbulence  statistics  in  these  simulations  compare 
very  well  with  DNS  results  of  Kim  et  al.  and  experimental  results  of  Nishino  and  Kasagi 
For  large  eddy  simulations  the  standard  Smagorinsky  model  with  constant  C*  =  0.1  and  van  Driest 
wall  damping  was  used  as  well  as  the  backscatter  model  (6)  with  two  test  cutoff  wavenumbers 
kc  =  l/2km  and  kc  =  SI 4km,  where  kc  =  akm  signifies  filtering  in  the  two  horizontal  directions 
with  cutoffs  kxc  =  akxm  and  kyc  =  akym,  respectively.  No  filtering  in  the  vertical  direction  was 
employed  in  the  models.  The  simulations  were  run  until  a  statistically  stationary  state  was  reached 
and  continued  for  additional  several  thousand  time  steps  to  collect  turbulence  statistics. 

Turbulent  rms  velocities  computed  for  all  models  are  plotted  in  Fig.  3  where  they  are  compared 
with  DNS  results  of  Gilbert  As  expected  from  numerous  tests  of  the  Smagorinsky  model 
reported  in  the  literature  it  predicts  rms  velocities  reasonably  well.  For  a  chosen  value  of  Gs  =  0.1 
the  streamwise  velocity  is  predicted  very  well  while  the  spanwise  and  normal  velocity  conaponents 
are  slightly  under  predicted.  The  net  effect  on  the  rms  velocities  of  adding  the  modelmg  expression 
(6)  to  the  Smagorinsky  model  is  small  since  the  models  were  designed  not  to  affect  significantly 
the  energetics  of  the  system.  The  backscatter  model  with  kc  =  3 /4km  shifts  peaks  of  the  horizontal 
velocity  components  too  close  to  the  wall  but  improves  the  normal  velocity  prediction  in  the  wall 
region.  For  kc  =  l/2km  the  streamwise  and  normal  velocities  are  similar  to  the  predictions  of  the 
Smagorinsky  model,  but  the  prediction  of  the  spanwise  component  is  significantly  improved  in  the 
wall  region.  This  case  provides  the  best  overall  improvement  in  the  rms  values.  The  backscatter 
model  (6)  has  also  been  implemented  using  a  top  hat  filter  with  the  filter  width  twice  the  mesh 
size  in  the  streamwise  and  spanwise  directions.  In  this  case  the  effects  of  the  backscatter  model 
were  found  to  be  very  small  and  the  results  for  the  rms  velocities  were  virtually  the  same  as  for 
the  pure  Smagorinsky  model.  The  differences  between  spectral  and  top  hat  filtering  applied  to 
the  backscatter  models  can  be  better  understood  by  calculating  a  nonconservative  part  of  the  SGS 
energy  flux  Fsgs  (or  SGS  dissipation  csas)  associated  with  the  models: 
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Figure  3:  Fluctuating  rms  velocities  normalized  by  Ut  as  a  function  of  the  distance  from  the  wall 
in  wall  units:  solid  line,  DNS  results;  broken  line,  Smagorinsky  model;  dotted  line,  Smagorinsky 
model  with  backscatter  model  (6)  for  kc  =  l/2km’,  broken-dotted  line,  Smagorinsky  model  with 
backscatter  model(6)  for  kc  =  SjAkm-  (a)  streamwise  component  Urma',  (b)  spanwise  component 
Vrma]  (c)  normal  component  Wrms- 
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FsGS  =  —^SGS  =  Tij 


dui 

'd^f 


(9) 


Contributions  to  this  quantity  from  the  Smagorinsky  model  and  the  backscatter  model  (6)  for  spec¬ 
tral  and  top  hat  filters  are  plotted  in  Fig.  4.  In  every  horizontal  plane  the  SGS  flux  is  decomposed 
into  its  negative  (forward  transfer)  and  positive  (backscatter)  components.  The  Smagorinsky  model 
is  purely  dissipative  providing  only  forward  transfer.  There  is  formally  no  distinction  between  re¬ 
solved  velocity  fields  for  both  filters  used  and  consequently  the  Smagorinsky  SGS  flux  is  the  same  in 
both  cases.  However,  different  test  filters  used  in  expression  (6)  will  give  different  results.  By  design 
the  backscatter  model  provides  approximately  equal  amounts  of  forward  transfer  and  backscatter, 
but  effects  of  the  model  calculated  with  the  top  hat  filter  are  significantly  smaller  than  for  the  sharp 
spectral  filter.  This  explains  its  small  effect  on  the  energetics  observed  in  the  rms  results.  The  sig¬ 
nificant  reduction  of  the  backscatter  for  the  top  hat  filter  as  compared  with  the  sharp  spectral  filter 
is  consistent  with  a  priori  analysis  of  Piomelli  et  al.  However,  since  the  resolved  scales  in  LES 
reported  in  this  paper  are  obtained  by  a  spectral  truncation  the  choice  of  the  spectral  filter  for  the 
backscatter  model  is  preferred. 


4  Conclusions 


We  have  used  the  properties  of  subgrid-scale  nonlinear  interactions  observed  in  DNS  of  turbulence 
to  propose  models  accounting  for  the  backscatter  phenomenon.  The  modeling  expressions  can  be 
added  to  any  standard  dissipative  SGS  model  improving  physical  realism  of  a  LES  procedure.  The 
models  provide  approximately  equal  amounts  of  forward-  and  back-scatter  and  affect  only  weakly 
the  energetics  of  the  system  as  predicted  by  the  dissipative  part  of  the  full  model.  The  models 
have  a  form  similar  to  the  nonlinear  term  in  the  Navier-Stokes  equations  and  thus  avoid  numerical 
stability  problems  associated  with  modeling  of  backscatter  by  a  negative  eddy  viscosity.  The  above 
properties  of  the  models  were  suggested  by  the  previous  a  priori  analysis  of  DNS  results  and  have 
been  confirmed  in  actual  large  eddy  simulations  for  turbulent  isotropic  and  channel  flows.  The  only 
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adjustable  paxameter  in  the  models  is  a  filter  width  which  provides  limited  control  over  the  amount 
of  backscatter. 
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